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Deconvolution  by  homomorphic  and 
Wieuer  Altering 

P.  Nicolas 


Executive  Summary:  The  Remand  for  more  sensitive, per<  cptton  of  sub¬ 
marine  signals  buried  within  ocean  noise  requires  statistical  methods  of 
analysis  A  statistical  sonar  theory  is  concerned  with  the  development  of 
probability  models  for  signals,  interferences,  and  underwater  experiment 
conditions,  and.  based  on  these  model", the  development  of  methods  for  the 
detection,  identification,  and  classification  of  submarines 

Studies  dealing  with  propagation  in  shallow  wa‘er  generally  model  the  re¬ 
ceived  signal  as  a  convolution  between  a  transmitted  pulse  (or  wavelet)  and 
the  medium  response.  In  this  focal  report  of  one  such  study  the  principal 
aim  is  to  extract  more  information  on  the  medium  -  such  as  backscattering 
effects  «uu  multipath  structure  uviTi  a  signal  received  cm  a  point  ■«<  rim  ui 
an  array  at  a  lower  signal-to  ncise  ratio  -  than  has  been  achieved  previously. 
This  clearly  could  have  a  direct  impact  on  future  sonar  system?. 

The  principal  advantage  of  the  so-called  method  is  that  it  does  not  require 
the  usual  Lssumption  of  minimum  phase  signals  (or  that  all  signals  have  a 
well  behaved  phase  structure)  and  is  il.cmuie  capable  of  coping  with  more 
realistic  propagation  conditions  where  in  general,  the  various  signal  arrivals 
have  a  complex  mixed-phase  structure. 

The  performance  of  the  method  is  demonstrated  using  both  simulated  and 
real  at-sea  data.  With  the  simulated  data,  deconvolution  of  the  wavelet  can 
b«  achieved  down  to  a  signal-to-n  oise  ratio  of  -5  dB,  while  the  n.ultipsths 
arc  well  separated  at  a  sigcai-to-noise  ratio  of  5  dB.  Using  an  explosive 
source  and  a  vertical  artsy  receiver  at  sea  one  can  separate  the  very  close 
reflected  and  refracted  paths  near  the  surface  in  the  order  of  1  or  2  ms. 
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Deconvolution  by  homomorphic  and 
Wiener  filtering 

P.  Nicolas 


Abstract:  This  study  is  concerned  with  deconvolution  methods  applied 
to  underwater  propagation  in  shallow  water,  whereby  the  received  signal  is 
modelled  as  the  convolution  between  the  transmitted  pulse  and  the  medium 
impulse  response.  The  aim  of  the  method  is  to  extract  information  on 
backscattering,  travel  time  delays,  boundary  reflection  and  refraction  from 
the  received  signal  on  a  point  receiver  or  an  array  for  both  seismic  and 
active  sonar  data.  Since  experimental  data  are  generally  mixed  phase,  due 
in  part  to  the  multiple  reflections  (bottom  and  surface),  the  conventional 
linear  filtering  which  assumes  the  minimum  phase  property,  loses  in  efficacy. 
In  order  to  handle  this  mixed  phase  characteristic  of  the  data,  we  proceed 
in  two  steps.  We  first  apply  a  homomorphic  filter  (complex  copstrum)  to 
deconvolve  the  wavelet.  Then  we  deconvolve  the  medium  impulse  response 
by  means  of  Wiener  filter.  The  efficacy  of  the  method  is  shown  on  both 
simulated  and  real  data  for  explosive  and  active  sonar  data,  t  •:  / 

.  ,  V-  •  . 

r  ) 
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1.  Introduction 

In  many  fields  of  physics,  such  as  geophysics,  seisinics  and  sonar,  we  are  faced  with 
pioblems  of  deconvolution.  The  observed  signal  received  from  a  sensor  in  these 
fields  is  often  considered  to  be  formed  by  the  convolution  of  the  transmitted  signal 
with  the  propagation  medium  impulse  response.  The  goal  of  deconvolution  filter¬ 
ing  is  to  recover  the  medium  impulse  response  from  the  recorded  signal.  Different 
methods  have  evolved  according  to  the  type  of  a  pnort  information  included  in  the 
signal  modelling.  If  the  transmitted  signal  is  known  exactly,  Wiener  filtering  is  con¬ 
ventionally  applied  under  the  assumption  of  minimum-phase  signals'  ( 1  —3) .  If  the 
source  signal  is  not  known  exactly  (which  is  the  case  for  an  explosive)  but  can  be 
modelled  by  a  parametric  transfer  function,  linear  prediction  methods  can  be  used 
with  success  [4—6]  However  these  methods  require  the  minimum  phase  condition. 
In  real  life  the  received  signal  is  generally  mixed-phase  which  is  the  case  for  seismic 
data.  When  considering  this  real  constraint,  another  approach  is  non-linear  filter¬ 
ing  based  on  the  generalized  superposition  principle  proposed  by  Oppenheim  and 
called  homomorphic  deconvolution  (7],  This  is  based  on  the  separation  of  the  so- 
called  wavelet  and  the  medium  impulse  response  in  the  cepstral  domain  [8,9].  Here 
we  present  a  method  which  combines  linear  and  non  linear  filtering  (10).  The  aim 
of  this  method  is  to  extract  information  on  backscattering,  travel-time  delays,  and 
boundary  reflection  and  refraction  from  the  received  signal  at  a  point  receiver  or  an 
array-  for  both  seismic  and  active  sonar  data.  Since  no  assumption  of  minimum 
phase  is  made,  we  first  apply  a  homomorphic  technique  (complex  cepstrum)  in  order 
to  deconvolve  the  wavelet.2  The  deconvolved  wavelet  is  then  taken  as  the  known 
signal,  and  we  estimate  the  boundary  reflections  and  travel  time-delays  by  means  of 
Wiener  filtering. 

The  report  is  structured  as  follows:  first,  we  define  the  wavelet  and  the  modelling 
of  the  medium  behaviour;  second,  we  advance  the  concept  of  homomorphic  de- 
convolution  and  its  mathematical  formulation;  third,  we  apply  Wiener  filtering  to 
the  recovery  of  boundary  reflection  and  propagat  ion  time-delays;  fourth  we  propose 
an  improvment  of  the  deconvolution  method  bared  on  a  combination  of  homomor¬ 
phic  deconvolution  and  Wiener  filtering;  and  fifth  the  application  to  seismic  and 
active  sonar  experiments  is  illustrated  We  present  results  obtained  on  both  simu¬ 
lated  and  field-recorded  marine  seismic  data  and  active  sonar  data.  We  point  out 
how  the  method  can  be  used  succcsfully  in  active  sonar  to  analyse  backscattering 
statistics  The  important  notion  of  minimum-phase  signals,  phase  unwrapping  and 
mathematical  investigation  of  the  complex  cepstrum  through  models  are  expanded 
in  appendices. 


The  term  ‘minimum-phase  signal'  is  defined  in  Appendix  A. 

The  term  wavelet  was  introduced  by  among  other  people,  Tribolet  [  1 1] ,  and  is  defined  and 
explained  in  the  first  part  of  the  present  report. 
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2.  Definition  of  the  wavelet 

This  section  is  devoted  to  the  definition  of  the  wavelet.  It  will  be  shown  later  that 
the  wavelet  contains  information  on  the  probability  characteristics  of  reverberation 
and  propagation  conditions.  The  study  of  the  statistical  features  of  reverberation 
presents  two  points  of  specific  interest:  one  is  the  properties  of  reverberation  as 
sonar  interference;  the  other  is  reverberation  as  a  phenomenon  which  helps  us  to 
estimate  the  properties  of  the  water  medium  and  its  boundaries.  Reverberation 
can  be  classified  into  three  types:  volume  reverberation,  reverberation  from  a  layer; 
and  reverberation  from  a  boundary.  A  propagation  signal  can  originate  from  an 
explosion  or  can  be  a  controlled  pulse  transmitted  from  a  point  source  and  received 
on  an  array  of  hydrophones  (here  a  vertical  ar-ay).  As  it  propagates  through  the 
medium  it  follows  three  paths:  the  direct  path;  thi  surface  reflected  path;  and  the 
bottom-layer  reflected  path.  Figure  1  represents  the  case  of  a  source  and  an  array 
closer  to  the  surface  as  opposed  than  in  Fig.  2,  which  shows  a  source  and  an  array 
closer  to  the  bottom.  These  figures  present  a  simplified  propagation  model  and  do 
not  take  into  account  the  ghosts  and  multiple  arrivals;  in  a  more  realistic  model 
these  can  be  removed  by  adaptive  linear  filtering  [12]. 

The  backscatterings  at  the  sea  surface  and  at  the  layer  boundary  are  defined  respec¬ 
tively  by  the  impulse  response  functions  and  MO-  The  medium  propagation 
is  defined  by  the  impulse  response  hm(t).  These  three  impulse  responses  are  ran¬ 
dom  processes.  In  the  first  case  (Fig.  1),  the  received  signal  is  dominated  by  the 
direct  arrival  and  the  surface-reflected  arrival.  The  layer-bottom  reflected  arrival 
comes  much  later  and  is  therefore  not  included  Assuming  that  the  medium  and  the 
surface  boundary  act  as  linear  filters,  the  signed  yi(f) — composed  of  the  direct  and 
surface- reflected  paths — is  given  in  the  time  interval  [0,T]  by 

yi(0  =  M(0  *  *(0  +  Am (0  *  MO  *  H*  -  T.)  *  *(0>  (1) 

where  z(t)  is  the  transmitted  signal  and  r,  the  propagation  time-delay  along  the 
surface  path.  In  the  second  case  (Fig.  2)  and  under  the  same  assumption,  the  signal 
yj(0 — composed  of  the  direct  and  bottom-layei  reflected  arrivals — is  given,  in  the 
time  interval  [O.Tj,  by 

y»(0  =  M(0  *  *(0  +  M(0  *  MO  *■  *(*  -  n)  *  *(0,  (2) 

where  rj  is  the  propagation  time-delay  along  the  bottom-layer  path. 

Because  the  impulse  responses  Am(0.M0  *nd  MO  random  processes,  yi(t) 
and  y*(0  defined  on  (0,7’)  are  to  be  considered  as  particular  realisations  of  random 


-  2  - 


SaCLAI  TCEN  SM-203 


signals.  By  taking  the  fourier  transform  of  both  sides  of  (1 )  and  (2),  we  have 

Yx(f)  =  Hm(f)X(f)  +  (3) 

Y,(f)  =-  HJf)X{f)  +  (4) 

where  Y\(f)  and  Vj(/)  are  particular  realiraticns  of  the  spectrum  of  the  signals  y,(t) 
and  yj(<)  respectively.  Both  of  the  equations  (3)  and  (4)  can  be  factored  in  two  ways 


Yy{f)  =  Hm(f)X(f)[  1  + 

] 

(5a) 

=  Hm(f)H.(f)X(f)e~>utT‘[  1  + 

H.(f) 

].  (6b) 

W)  =  f/m(/)X(/)[l  +  f/,(/)e-^T' 

1 

(6a) 

=  Hm(f)Hi{f)X{f)e-^tT'[  1  + 

_ 1 _ 

(6b) 

Hi(f) 

By  definition,  the  wavelet  is  given,  with  respect  to  the  factored  expression  by 


Hm(f)X(f)  (first  factored  expression) 

(second  factored  expression), 


VK, 


(first  factored  expression) 
(second  factored  expression). 


Depending  on  the  factored  expression,  the  wavelet  contains  information  on  only 
medium  propagation  or  on  both  medium  propagation  and  boundary  backscattering. 
In  order  to  separate  the  wavelet  from  the  other  components,  we  can  take  the  complex 
logarithm  of  Y\(f)  and  yj(/).  If  the  surface  impulse  response  is  minimum  phase 
(definition  and  details  on  minimum-phase  signals  are  given  in  Appendix  A),  the 
modulus  of  H,{f)  is  less  than  unity  and  one  uses  the  first  factored  expression. 

logy,)/)  =  logtfm(/)X(/)[l  +  H.(f)e-^'r>] 

=  log  DM/)  +  tf*(/)e->2w/T‘  -  U.(/)5e-'WT*. 

If  the  surface  impulse  response  is  not  minimum  phase,  the  module  of  H,(f)  is  greater 
than  unity  and  the  second  factored  expression  is  used: 

log  i'll/)  =  logH„(/)log«,(/)A:(/)«-',-'r'|l  + 

=  loglV,(/)+  -o'”"'. 
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We  can  do  the  same  for  the  bottom-layer  reflected  path.  Then,  by  using  an  appro¬ 
priate  Alter  and  suppressing  the  linear  phase  component  when  the  impulse  response 
is  not  minimum  phase,  we  can  extract  the  wavelet  from  the  received  signal.  In  fact, 
instead  of  Altering  the  signal  in  the  frequency  domain,  we  Alter  the  inverse  fourier 
transform  of  the  complex  logarithm,  which  is  the  complex  cepstrum  by  deAnition 
(see  Sect.  3). 

In  brief,  the  wavelet  is  an  artiAci&l  transmitted  signal  in  the  sense  that  it  represents 
the  transmitted  signal  modified  by  the  propagation  and  backscattering  characteris¬ 
tics.  Depending  on  the  boundary  properties,  the  wavelet  carries  more  information  or 
less  information  (minimum  or  not  minimum  phase  property).  By  extracting  the  im¬ 
pulse  response  functions  hm(t),  h,(t)  and  h|(<),  we  improve  the  modelling  of  medium 
propagation  and  backscattering  statistics  which  can  be  compared  to  existing  theo¬ 
retical  models. 
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3.  Homomorphic  deconvolution 


3.1.  HOMOMORPHIC  SYSTEMS 

In  reverberation  we  are  faced  with  the  problem  of  filtering  signals  that  have  been 
combined  by  convolution.  It  would  be  avantageous  to  transform  these  non-linear 
systems  into  linear  systems  by  applying  the  appropriate  filtering.  This  leads  to 
systems  which  obey  the  ‘generalized  principle  of  superposition’.  Given  two  inputs 
let  us  assume  that  they  are  related  together  by  a  rule  o.  If  a  is  a  scalar  let  :  be 
a  rule  to  combine  s  with  any  of  the  two  inputs.  Similarly,  we  denote  o  the  rule  to 
combine  the  outputs  together  and  •  a  rule  to  combine  a  scalar  with  an  output.  If 
H  is  the  system  transformation,  we  state: 

:*,(«))  = 


The  systems  that  verify  the  two  preceding  equations  are  said  to  obey  a  ‘generalized 
principle  of  superposition’  [9].  It  can  be  shown  that  if  ihe  system  inputs  constitute  a 
vector  space  with  the  operations  o  and  :  corresponding  to  vector  addition  and  scalar 
multiplication  and  the  system  outputs  constitute  a  vector  space  with  the  operation  o 
and  •  corresponding  to  vector  addition  and  scalar  multiplication,  then  all  systems  of 
this  kind  can  be  represented  as  a  cascade  of  three  systems  referred  as  the  ‘canonical 
representation  of  homomorphic  systems’,  shown  in  Fig.  3. 

The  first  system  29  0  has  the  following  property: 

T>e[*i (t)  o  *2(0]  =  A>(*i(t)]  +  Do[*i{t)}  =  *i(0  +  *2 (0- 
D0[s  :  *i(t)]  =  jD0[*i(fj]  =  •’ii(t). 

The  effect  of  the  system  D„  is  to  transform  the  signals  ii(t)  and  £2(1)  according  to 
the  rule  o  into  a  conventional  linear  combination  of  corresponding  signals  Z70[xi(l)] 
and  D(,[xi(t)].  The  system  L  is  a  linear  system: 

L[xx(t)  t-  x2(0]  =  L[jt(0]  +  £[*2(0]  =  3/1  (0  +  jMO. 

I[sii(t)]  =  aX.(ii(t)]  =  ayt(t). 

The  system  D~ 1  transforms  from  addition  to  the  rule  o: 

D7l[y\(t)  +  y*(0)  =  Z?0-1(i/1(t)]oZ?0-l[y2(t)]  =  yx  (t)oyj(t), 

^oM^yOO]  =  a*  Dcl{yi(t)}  =  3  •  yx(t). 
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All  the  homomorphic  system?  with  the  same  input  and  the  sam.  output  differ  only 
in  the  linear  part.  Consequently,  by  choosing  the  transformations  D0  and  D0  we 
are  left  with  a  linear  problem. 

We  are  going  to  apply  these  results  to  convolved  inputs  signals.  The  rule  o  becomes 
the  convolution  *  We  choose  the  output  rule  o  to  be  equal  to  the  input  rule  and 
therefore  o  is  also  equal  to  tl  e  convolution  *.  The  canonical  representation  of  an 
homomorphic  deconvolution  system  is  shown  in  Fig.  4. 


3.2.  COMPLEX  OEPSTKUM 

■  3.2.1.  Mathematical  representation  of  the  system  D.  and  definition  of  the  complex 
cepscrurn 

The  system  D.  is  defined  by  the  property  that  the  2-transform  (or  the  fourier  trans¬ 
form  on  the  unit  circle)  of  its  output  is  equal  to  the  complex  logarithm  of  the 
2-transform  (or  the  fourier  transform  on  the  unit  circle)  of  its  input: 

i(n)  =  £>,|x(n)], 

X(z)  =  logX(j), 

where  x(n)  is  the  nth  sample  of  x(t).  According  to  this  definition,  the  characteristic 
system  D.  of  the  homomorphic  deconvolution  is  as  shown  in  Fig.  5. 

The  output  of  the  system  Dm,  denoted  i(n),  is  called  the  complex  cepstrum  of  the 
input  signal  x(n).  This  terminology  is  used  by  analogy  to  the  power  cepstrum  defined 
by  Bogert,  Healy  and  Tukey.  Specifically,  the  cepstrum  of  a  signal  was  defined  as 
the  power  spectrum  of  the  logarithm  of  the  power  spectrum. 


Remark  These  quantities  are  not  too  far  from  each  other,  because  the  power 
cepstrum  is  proportional  to  the  even  part  of  the  complex  cepstrum. 

■  3.2.2  Definition  of  the  complex  logarithm 

In  this  section  the  complex  logarithm  chosen  as  the  homomorphic  system  D.  is 
defined.  One  first  sets  the  definition  of  the  logaritlim  and  then  considers  more 
particularly  the  phase  unwrapping  problem.  Its  prevalent  lole  and  the  critical  points 
of  the  different  phase  unwrapping  methods  are  pointed  out.  To  finish,  some  examples 
of  phase  unwrapping  are  given. 


■  3.2  2a.  Definition 

Let  be  x(n)  a  real  sequence  and  X(z)  its  r-tianaform.  one  wants  to  define  the 
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logarithm  of  _Y(2).  The  complex  logarithm  is  a  ‘multi-valued  function’  and  therefore 
one  must  choose  a  determination  for  which  the  logarithm  is  a  continuous  function. 
Usually  one  takes  the  ‘principal  value  determination’  defined  by 

log-Y(z)  =  log  |JY(x)|  +  j  Arg[jY(z)), 

where  Arg(X(z)j  €  (~jt,  +  jt). 

All  the  otner  determinations  are  obtained  by  adding  a  multiple  of  2ir  to  Arg[X(2)| 
In  our  case  the  sequence  z(n)  is  a  convolution  of  two  sequences,  ati(n)  and  i2(n): 

logA(2)  =  logXi(z)  +  logJY2(z). 

The  principal  value  of  the  logarithm  oi  me  pioduct  of  complex  sequences  i9  not 
always  the  sum  of  the  principal  value  of  each  of  the  signals.  This  is  in  contradic¬ 
tion  with  the  unicity  of  the  homomorphic  system  D.,  and  means  that  the  complex 
logarithm  cannot  be  defined  from  the  principal  determination  alone.  Besides  given 
properties  of  the  sequences  z(n)  one  needs  another  definition  of  the  complex  loga¬ 
rithm. 

The  complex  logarithm  will  be  defined  from  its  derivative.  If  one  assumes  a  single¬ 
value  differentiable  complex  logarithm  (principal  value)  and  the  analyticity  of  X(z) 
one  can  derive  the  phase  as  follows: 

The  evaluation  of  the  complex  logarithm  on  the  unit  circle  z  =  e;u'  is  performed  in 
the  following  manner: 

_doM  =  1  dA(e^)  =  1  dA(e^)  1 

dz  '  j z  du>  JY(eJW)  du>  jz  ’ 
so 

dA(e^)  _  _ 1  dA(e^) 

dw  X(ei“)  dw 

Given  JY(e-,“)  —  A'/i(e;w)  +  jX/(e^u),  we  have 

dA(e>)  =  dXn(e^)  ^  .  dXf(e^) 

&u)  dw  dw  * 

x\^)  =  xR{*n  +  ix'Mn, 

where  the  prime  indicates  the  differentiation  with  respect  to  w. 
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Hence, 


jfV)  _  PW")  -f  j^PgK)  -  jXfen I 

X(ei“)  |AT(e>-)|2 

=  f  x't[enxi(*n\ 

h  j(X;(e^).Y«(e^)  -  X;(e-'-)X/(ei“)|] . 


we  have 


*V")  -  J(£)  =  ^loe|;c(c",i  + 


da,  gl  {  |JY(e>‘)|* 


Property  of  the  phase  derivative:  The  phase  derivative  is  an  even  function  of 
u  since  z(t )  is  a  real  function. 


Proof: 


and  hence, 


*(0  =  f  +  °°  x(t)e-^dt  =  XR(eiu)  +  jX[(eju) 

J  —  no 

X{e~jM)=  (  x(tyutdt  =  A«(e^)  4 -  j*/(e>“) 

v  —  OO 


=  xRy“)  -  jXtien 


X(e~^)=  XR(e-^)  +  jX/(e-^), 


X*(e-'")  =  A^e'""),  X,(e~^)  =  -Xj(e^), 

x'ni'-n  -  -Xr(^Y  x'fa-i")  - 


-^arg[A(e^))  =  arg[X(e^)j. 


Assumption:  both  X(z)  and  A(r)  are  analytic  in  a  region  included  the  unit  circle 
X{z)  and  A(z)  have  no  singularities  on  the  unit  circle.  Consequently  the  functions 
X/fe^"),  XR(e^w),  Xr(e-’u)  and  XR(e^w)  are  analytic  and  the  phase  derivative  is 
analytic  in  the  convergence  domain  of  X(z)  and  A(z).  Let  us  recall  the  following 
theorem: 
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Theorem:  Let  fl  be  8  continuum  of  the  complex  plane  end  /  a  continuous  function 
on  fi.  Tne  necessary  and  sufficient  condition  for  the  function  /  to  have  a  primitive 
is  that  the  integral 

Lf{z)dz 

is  null  for  any  contour  C  included  in  fi.  Under  this  condition  all  the  primitives  F(z) 
are  obtained  by  the  formula 

F(z)  =  I  /(u)du  +  K, 

•'*0 

where  is  a  point  of  fi  and  K  and  arbitrary  complex  constant.  /( u)du  is  the 
integral  of  /  on  any  path  of  fi,  starting  from  the  point  zq  and  joining  the  point  z. 

As  d/dz.Y(z)  is  analytic  in  the  convergence  domain  of  ^f(z)  and  A'(z),  we  have 
according  to  the  Cauchy  theorem, 

J  -~arg(A'(z)Jdz  =  0 

on  any  contour  C  included  in  the  convergence  domain.  According  to  previous  results 
the  phase  is  defined  without  any  ambiguity  on  the  unit  circle  to  within  an  integration 
constant: 

«g(  Ww)|  =  J"  arg[JT(e>")I  du>  +  K. 


The  constant  K  is  evaluated  in  the  following  way:  the  complex  logarithm  must  fulfill 
the  requirement,  given  two  sequences  zy  (n)  and  z2(n), 

log  Xi(e^)XJ(e}U)  =  logXx(e^)  +  log  Jf,(e>w) 


which  is  equivalent  to 

log  |Jrl(e*-)*,(e'“)|  +  j  arg[A'l(e^)A'2(e^)] 

=  log|JfI(^“)|  +  j  arg[Xi(e^)j  +  log)*,^)!  +  ;  arg[Jf,(e*-)]. 


One  must  have  arg(A'l(e;“')A2(eJU')]  =  arg[A'i(e''u')]  +  arg[A2(eJ“')),  or 

j  ^arg(A1(e^)A:(ejw)]du;  +  A'u  =  ^“^arg(A-,(e^))du;  +  A'1 

+  *rg[XJ(eiu)]dw  +  Ki. 

Jo 


To  have  this  equality  verified  for  any  sequences  *i(n)  and  r2(n),  the  constants 
Kt ,  Kj  and  K\ j  must  vanish.  To  have  the  constant  K  =  0  means  arg[J!f(e'7‘*')|w=0]  = 

0;  tut  J¥(ef-)|w=o  =  Ln=-oc  z(n)-  40(1  80  "glE^T-oo  z(n))  -  0 
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■  3.2.2b.  Properties  of  the  phase  arg[Jf(e-'“))  and  requirements  for  the  signal  z(n) 
The  determination  of  the  constant  K  leads  to  a  specific  property  for  the  sequence 
*(n):  the  dc  component  (polarity)  must  be  positive. 

And  the  phase  has  the  following  properties: 

(1)  arg(A(e>w)|  is  an  odd  function  of  u>  such  that 

arg[*(e-m  =  j  "  arg(X(e-^)j  du, 

=  j  -^arg[A(e-,u,)]du> 

-  -  arg[X(e-’")), 

(2)  arg(A' (e-,w)]  is  a  continuous  function  of  ui, 

(3)  arg[A (eJa')|u,=  IT]  -  /„”  d/dw  a^g(A'(e■''*,)]  du>  =  0  because  arg(A(e-’“)]  is  peri¬ 
odic  in  u)  with  a  period  2ir  and  is  an  odd  function  of  u>. 

Since  the  phase  derivative  is  an  even  function  of  w,  we  have 

/  •£;  arg[^(eJU,)]dw  =  -  y  arg[A(e"-?u’)J  du> 

=  l  arg(Jf  (c>“)]  dw. 

And  so, 

f  -Aarg*(e^)dw  =  -  f  -1  arg  X(e>“)  du>, 

2  JT  J  dw  7T  J0  dw 

and  the  previous  requirement,  for  u>  =  ir,  leads  to  a  second  property  of  the  sequence 
z(n):  z(n)  must  have  a  zero-mean  phase  derivative. 

Conclusion:  On  the  space  of  the  sequences  x(n)  with  a  positive  dc  component  (po¬ 
larity)  and  zero-mean  phase  derivative,  the  complex  logarithm  defines  an  invertible 
homomorphic  system. 

■  3.2.3.  Phase  unwrapping 


l  3.2.3 a.  Principle 

The  phase  unwrapping  involves  computing  a  continuous  phase  from  the  set  of  princi¬ 
pal  phase  samples.  Various  techniques  have  been  developed.  A  basic  one  is  Schafer’s 
Algorithm,  which  is  based  on  the  the  following: 
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(1)  whenever  a  jump  -2k  is  detected  while  unwrapping  along  the  positive  w-axis 
a  constant  2k  is  added  to  the  principal  value  at  that  point,  and 

(2)  whenever  a  jump  of  2k  is  detected  while  unwrapping  along  the  positive  w-axis 
a  constant  -2*  is  added  to  the  principal  value  at  that  point,  with  a  ’jump’ 
defined  as  the  difference  between  a  new  principal  value  and  an  old  one. 


A  jump  has  a  threshold  defined  in  terms  of  the  difference  between  the  two  principal 
values  at  adjacent  frequencies;  below  this  threshold  the  jump  does  not  exist.  This 
implies  a  frequency  sampling  fine  enough  to  set  the  difference  in  the  principal  values 
of  two  adjacent  samples  be  detected  as  a  jump.  However  although  this  algorithm  is 
simple  to  implement  it  does  not  provide  accurate  results  in  the  case  of  a  sharp  phase 
curve,  since  there  is  only  principal  value  phase  information  and  this  is  not  sufficient. 

To  overcome  this,  we  have  completed  the  phase  unwrapping  by  a  modified  Tribolet 
Algorithm,  which  takes  into  account  the  information  in  the  first  derivative  of  the 
phase.  Let  us  recall  briefly  the  principle  of  the  Tribolet  Algorithm  [11]. 


One  calculates  the  phase  at  the  frequency  w  from  the  mean  of  the  integral 


-)  -  f 


arg  X(en=  I  ^  arg  *(eJW)  dw. 


This  integral  is  approximated  by  the  trapezoidal  rule.  Assuming  that  the  phase  is 
known  at  the  frequency  w<,  one  estimates  the  phase  at  the  frequency  w1(  i  (w1+i  >  u<) 

by 


arg(A'(eJW'  +  l  )/w<]  =  arg  Xie^’)  4- 


aigX[e^‘*1)  +  -j^arg  X(eju>‘) 


where  Aw  =  Wj+ 1  —  Wj . 


A  phase  estimate  is  called  consistent  if 


3k(wi+})  B  |arg(X{ejw^‘)/w<j  -  Arg[Jf («'"*)]  +  2x*(w<+1)|  <  THLDl  <  #. 


The  idea  of  the  algorithm  is  to  adapt  the  step  size  Aw  until  the  phase  estimate  be¬ 
comes  consistent.  The  algorithm  requires  a  second  threshold  THLD2  in  order  to  con¬ 
trol  the  phase  increase  between  two  consecutive  frequency  samples  The  unwrapped 
phase  argJf(c,“’'+l)  at  frequency  /  is  used  to  estimate  the  phase  at  frequency  Wj+j 
and  so  on.  One  recalls  that  the  phase  derivative  is  given  by 

i  ,L1J  ... 
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Thus  it  can  be  computed  very  fast  using  the  FFT  according  that 

=  -j'FFT(n*(n)|. 


This  algorithm  works  well  as  long  as  the  spectrum  does  not  have  any  zeros  close  to 
the  unit  circle.  In  this  case  the  phase  derivative  given  by  the  previous  relationship 
and  computed  by  FFT  has  singularities  and  presents  big  spikes.  Thus  the  phase 
increase  is  no  larger  controlled.  To  improve  the  T^ibolet  algorithm,  an  idea  has 
been  suggested  by  (13).  It  consists  of  fitting  a  curve  to  the  phase  derivatives  before 
performing  the  numerical  integration:  one  fits  cubic  splines  S(u>),  having  continuous 
first  and  second  derivatives,  to  the  phase  derivative.  The  phase  is  then  given  by 

arg  X(e*u)  ~  /"  S(e'“)  dw, 

Jo 


or,  according  to  [14] 


arg(Jf(e-'“,+  1  )/<*>,)  -  aigX(e-,“' )  + 


r-^-argA'(e^+‘)+  -iarg^(e^) 
aw  aw 


-r-Sit***")-  —S(e’u') 


du> 


du> 


where 

This  can  be  also  computed  by  FF 1  as 

-  XI(ejw)XR(eju'))  +  2XR(enX,(en 

K^/V")  -  Xi(e>“))] 


and 

XR(^)  +  jXj{en  =  -FFT[n3*(n)]. 
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■  3.2.3b.  Examples  of  phase  unwrapping 
(a)  First  example 

Let  y(n)  be  a  time  series  which  is  the  convolution  of  two  time  series:  w(n)  which  is 
a  CW  pulse  windowed  by  a  Hanning  window  and  r(n)  given  by 

r(n)  =  <5(n  -  )  +  6(n  -  n j)  +  £(n  -  n3). 

The  following  processing  is  applied  to  the  time  series  y(n): 

(1)  compute  the  spectrum  of  y(n)  by  FFT; 

(2)  band-pass  filter  the  spectrum  around  the  CW  pulse  frequency; 

(3)  apply  the  band-pass  mapping  system; 

(4)  compute  the  first  and  second  derivatives  of  the  shifted  and  stretched  spec¬ 
trum; 

(B)  compute  the  unwrapped  phase. 

y(n)  is  a  256  time-sample  series.  The  normalized  CW  frequency  is  0.25  and  the 
Hanning  window  length  is  64  time-samples.  The  time  series  r(n)  is  given  by 

r(n)  =  *(n  -  55)  +  -  90)  +  6(n  -  125). 

The  band-pass  mapping  is  defined  as  in  iubsubsect.  3.2.6: 

!y|e(->“J))  ^  o,  for  uq  <  |u>|  <  u>2, 

F*[e<^>]^0,  for  M=u,1(u;2, 

0,  otherwise. 

The  spectra  of  tu(n),r(n)  and  y(n)  after  band-pass  mapping  are  depicted  in  Figs.  9, 
10  and  11  respectively.  The  phases  of  u/(n),r(n)  and  y(n),  before  unwrapping,  are 
represented  in  Figs.  12,  13  and  14.  The  first  and  second  derivatives  of  the  CW  pulse 
phase  are  represented  in  Figs.  15  and  16.  The  first  and  second  derivatives  of  the 
medium  response  phase  are  represented  in  Figs.  17  and  18.  The  first  and  second 
derivatives  of  the  received  signal  phase  are  depicted  in  Figs.  19  and  20  respectivelv 
The  unwrapped  phase  of  the  wavelet,  the  medium  response,  and  the  received  signal 
before  removed  of  the  Linear  phase,  are  repiesented  in  Figs.  21,  22  and  23  respectively. 
The  unwrapped  phase  of  the  wavelet,  the  medium  response,  and  the  received  signal 
after  removal  of  the  linear  phase,  are  represented  in  Figs.  24,  25  and  26  respectively. 


Remark  The  band-pass  mapping  introduces  some  small  instabilities  into  the  phase 
around  the  cut-off  frequencies  uq  and  u>2 .  The  instabilities  are  well  shown  on  the 
first  and  second  phase  derivatives  of  the  wavelet  phase. 
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(b)  Second  example 

Field  marine  explosive  data  have  been  recorded  at  the  output  of  a  vertical  array  of 
32  hydrophones.  The  received  signal  path  is  composed  of  a  direct  path  followed  by 
a  refracted  and  a  reflected  path  at  the  sea  surface.  The  bottom  interaction  comes 
much  later  and  is  not  accounted  for  in  the  present  data.  One  is  looking  at  the  output 
of  thrt  hydrophone  17.  The  phase  unwrapping  is  processed  on  the  full  frequency  ?iand 
(no  band-pass  mapping)  and  the  results  are  shown  in  Figs.  27  and  28. 


■  3.2.4.  Properties  of  the  complex  cepstrum 

The  complex  cepstrum  has  some  properties  which  are  useful  for  the  design  of  filters 
and  transmitted  signals.  Some  of  these  properties  are  summarized  below. 


Property  I  The  complex  cepstrum  of  a  convolution  of  two  (or  more)  signals  is 
the  sum  of  the  individual  complex  cepstra. 

Property  2:  The  complex  cepstrum  y(n)  of  a  minimum  phase  sequence  y(n)  is 
zero  for  n  <  0,  and  the  complex  cepstrum  of  a  maximum  phase 
sequence  is  zero  for  n  >  0.  (See  the  definition  of  a  minimum  and 
maximum  phase  sequence  in  Appendix  A). 

Property  3:  The  complex  cepstrum  of  a  pulse  whose  spectrum  is  smooth  tends 
to  be  concentrated  around  low  frequency  values. 


Property  4:  The  complex  cepstrum  of  a  periodic  impulse  train  is  a  periodic 
impulse  train  with  the  same  period. 


■  3.2.5.  Sensitivity  of  the  compiex  cepstrum  to  the  noise 

The  more  critical  part  of  the  complex  cepstrum  is  the  unwrapping  of  the  phase  due 
to  its  sensitivity  to  the  additive  noise.  In  the  following  discussion  we  try  to  show 
how  the  behaviour  of  the  signal  phase  depends  of  the  signal- to- noise  ratio  and  the 
noise  phase.  The  received  signal  plus  additive  noise  can  be  expressed  as 

a(t)  =  y(f)  +  n[t), 

where  y(t)  is  the  convolution  of  two  or  more  signals  and  n(t)  is  the  additive  noise. 
In  the  frequency  domain  this  equation  becomes 

S(u>)  =  Y(u>)+N(v) 


and 

log5(w)  =  Iog[r(«)  +  JV(a/)]. 


Equation  (7)  can  be  rewritten  as 


logS(w)  =  logF(w)  +  log 


W 


(7) 
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If  we  consider  a  signal-to-noise  ratio  which  is  relatively  high,  we  can  assume  that 

Y(u)  <  l' 

and  Eq.  (7)  can  be  expanded  into  its  Taylor  series  as  follows: 

/V(w)  l/tf(c)\S 


logS(u>)  =  log  Y  (a;)  + 


Y(u>) 


_  i  wy 

2  Wm; 


+ . . .. 


If  we  consider  the  terms  of  the  series  of  the  second  order  and  higher  to  be  negligible 
the  phase  of  S(u>)  can  be  expressed  as  follows: 


*y(w)  +  sin[$N(n>)  - 


(8) 


where  $n(w)  and  $y(“^)  we  respectively  the  phase  of  the  noise  and  the  signal.  If 
we  now  consider  that  the  signal-to-noise  ratio  is  low  such  that 


A^) 

>» 


> 


1, 


Equation  (7)  can  be  approximated  as 


log5(u/)  =  log  jV(u>)  + 


W 


and  the  phase  of  S(w)  is  equal  to 


*h(«)  +  ^  sin[§YM  ~  $n(w)1.  (9) 

Equations  (8)  and  (9)  show  that  the  phase  of  the  received  signal  s(t)  can  become 
unpredictable  (random)  because  of  the  noise.  When  the  signal-to-noise  ratio  is  low 
(Eq.  9)  the  phase  is  dominated  by  the  phase  of  the  noise.  The  part  of  the  noise 
spectrum  which  is  not  overlapped  by  the  signal  spectrum  can  be  removed  by  band¬ 
pass  filtering  in  order  to  avoid  the  situation  of  a  low  signal-to-noise  ratio. 


Remark  The  use  of  band- pass  filters  before  the  homomorphic  deconvolution  leads 
to  the  notion  of  band-pas3  systems  (see  Subsubsect.  3.2.6). 
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■  3.2.6.  Definition  of  the  band-past  mapping  system 

In  many  applications  the  signals  have  band- pass  characteristics.  In  general  the 
signals  are  band-pass  filtered  before  being  sampled  in  order  to  increase  the  signal-to- 
noise  ratio.  The  homomorphic  deconvolution  as  described  above  cannot  be  applied 
directly  to  the  band-pass  signal:  the  logarithm  is  not  defined  in  the  frequency  domain 
where  the  signal  vanishes.  Before  applying  any  cepstrum  analysis  one  must  find  a 
system  which  transforms  the  band-pass  signal  into  a  full-band  signal.  Such  a  system 
is  called  a  band-pass  mapping  system. 


a  3.2.6a.  Principle 

This  notion  of  a  band-pass  mapping  system  has  been  introduced  by  TVibolet  [11]. 
Let  *(n)  be  a  stable  sequence  and  X (e^l‘d)  its  fourier  transform  satisfying 


!  Y(e^“d)  ^o,  for  <  |w|  <  u»j 
**[e(>'")|  /  0,  for  M  =  wi.uij 
0,  otherwise, 

where  u/j.wj  are  the  cut-off  frequencies.  Let  BP  denote  the  band-pass  mapping 
system  operator  defined  by 

i(n)  =  BP(z(n)] 

such  that  the  fourier  transform  of  i(n)  verifies 


where 


X[e(j,:;):  =  *[e<'w>),  0  <  ja>i  <  t, 


u>(u;)  =  t 


U>  —  U>1 

-  uj\  ’ 


£  |u/|  S 


(10) 


Remark  This  frequency  transformation  is  a  scaling  operation  that  shifts  and 
stretches  the  signal’s  pass-band  to  occupy  all  of  the  frequency  band. 


TVibolet  has  verified  that  this  band-pass  mapping  is  an  invertible  homomorphic 
operation  with  convolution  as  input  and  output  operations.  The  inverse  operation 
is  defined  by 

i  _  Ulj  —  ta)j  i 

ui  =  4-  ui - ,  wj  <  <  u>j, 


for  <  |u|  <  u>j, 
otherwise. 


From  Eq.  (10),  we  have 


=  X[e(>w)j,  forw  =  u/, 
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and  thus 

x'(n)  =  z(n). 


The  block  diagram  of  the  band-pass  complex  cepstrum  system  D.  is  illustrated  in 
Fig.  6.  The  band-pass  mapping  system  is  illustrated  in  Fig.  7  using  simulated  data. 
The  received  signal  is  the  convolution  of  a  Hanning  windowed  CW  pulse  with  three 
Dirac. 


■  3.2.6b.  Implementation  of  the  band-pass  mapping  system 

Let  x(n)  be  a  N  (power  of  2)  samples  sequence  and  X(n)  the  corresponding  DFT 
sequence  One  may  suppose  the  sampled  spectrum  to  be  symmetrically  band-pass 
filtered  around  the  normalized  frequency  0.5,  which  corresponds  to  the  frequency 
sample  J-A;  the  sampled  cut-off  frequencies  are  ^ N-N j  and  J-JV  4-  .The  band-pass 
mapping  operation  can  be  decomposed  into  the  four  following  phases: 

(1)  shift  the  band-pass  spectrum  to  0, 

(2)  compute  the  2 Nx  +  1  inverse  DFT  of  the  sequence  X(n)  for  f  A  -  TVj  <  n  < 
£  4-  Ai  in  order  to  get  a  2Ai  +1  time  s-'-ies, 

(3)  zero-pad  this  new  time  series  to  get  a  time  sequence  of  N  samples, 

(4)  compute  the  N  inverse  FFT. 

Operations  2,  3  and  4  correspond  to  the  stretching  of  the  spectrum  (spectrum  in¬ 
terpolation). 

The  inverse  band-pass  mapping  operation  can  be  broken  down  into  the  four  following 
phases: 

(1)  cut  the  N  samples  deconvolved  sequence  at  the  first  2At  \  samples; 

(2)  compute  the  2 Nx  +  1DFT  of  the  2Ai  +  1  sequence; 

(3)  shift  the  spectrum  to  the  ;j;A  -  Nx  frequency  sample;  set  the  spectrum  value 
X(n)  at  0  for  1  <  n  <  J-A  -  Nx  -  1  and  -f  JV,  +  1  <  n  <  f  A ; 

(4)  compute  the  N  inverse  FFT. 


■  3.2.7.  Normalization  of  the  signal  before  applying  the  complex  cepstrum 


m  3  2.7a.  Principle 

When  the  time  series  i(n)  does  not  fulfill  the  requirements  that  its  dc  component 
(polarity)  is  positive  and  its  mean  phase  derivative  is  equal  to  zero,  the  input  se¬ 
quence  £(n)  must  be  normalized  in  order  to  be  able  to  app.y  the  complex  cepstrum. 
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We  recall  that  if  *(n)  has  a  mean  phase  derivative  no  null, 

-  /  -^-argX(e°w)]dw  =  -  arg  X ie'Jw)]|u,=, 
r  Jo  aw  x 

The  first  part  of  the  normalization  consists  of  multiplying  X[etJU,l]  by  e* -J"r\ 


=  r. 


I  [*’  d 

x  J0  du- 


arg 


*[e( 


->"r) 


dw  =  -argX(e°w))  -  -  / 

X  X  J  Q 


r  Qwi  =  t  -  r  =  0. 


The  second  part  of  the  normalization  consists  of  multiplying  X [e^"1]  by  the  polarity 


■  3.2.7b  Restoration  of  the  linear  phase  components 

Consider  a  signal  x(t)  which  is  the  convolution  of  two  signals  *i(t)  and  xj(t)  and 
the  respective  spectra  for  which  are  given  by 

X,(f)  -  A'^n[(/)^WT,  , 


where  X\ nl(f)  and  Xi„i{f)  are  respectively  the  non-linear  phase  components  of 
Xi  (/)  and  Xj (/).  Then  take  the  logarithm  of  X(f) 

logX(f)  =  log  |X, „,(/)|  +  j  arg  X\„i(f)  +  log  |*Jr,((/)l 
+  j  arg Xj„i{f)  ~  j2xfr\  -  j2xfr2. 

And  then  remove  the  the  linear  phase  component: 

JogA„((/)  -  log •  Alnjf /)|  +  j  argX)ni(f)  +  log |.Yln/(/)|  +  j  arg XJni(f). 

This  last  relation  shows  that  we  get  an  infinity  of  solutions  comprising  all  the  signals 
with  the  same  non-linear  phase  components.  Assuming  that  we  are  able  to  separate 
ilni  and  *jnj  in  the  cepstral  domain,  we  must  restore  the  proper  linear  phase  to 
each  of  the  deconvolved  signals  iIr,i  and  ij „<•  This  task  becomes  infeasible  if  we  do 
not  have  a  prion  information  on  the  original  signals  xt(t)  and  x2(t).  For  example, 
if  we  assume  that  one  of  the  signals  has  no  linear  phase  component,  let  us  say  xj(t), 
then  rj  equals  zero  and  the  global  linear  phase  is  restored  to  it(t).  At  this  point  it 
is  more  a  matter  of  experimental  conditions,  as  we  can  see  in  the  example  treated  in 
Appendix  B.  In  the  results  concerning  the  active  sonar  simulation,  the  deconvolved 
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wavelet  is  rescaled  in  time  by  computing  the  cross-correlation  with  the  transmitted 
pulse. 

The  same  problem  arises  for  the  signal  polarity. 

The  global  complex  cepstrum  deconvolution  system  is  depicted  in  Fig.  8a. 
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4.  Deconvolution  of  the  medium  response  by  Wiener  filtering 

The  goal  of  this  section  is  to  provide  a  method  of  deconvolving  the  medium  impulse 
response.  Since  the  homomorphic  deconvolution  was  not  the  best  one  for  estimating 
the  medium  response  (see  Appendix  B)  in  the  presence  of  additive  noise,  we  use 
a  digital  Wiener  filter  It  belongs  to  the  class  of  linear  time-invariant  filters  with 
a  criterion  of  minimization  of  the  mean  quadratic  error  [27,28].  It  attempts  to 
optimally  transform  a  given  signal  to  another,  here  the  received  signal,  into  the 
medium  impulse  response.  It  is  stressed  that  the  estimation  of  the  length  and 
the  lag  of  the  Wiener  filter  are  not  treated  here;  those  are  discussed  in  [16].  The 
Wiener  filter  assumes  that  we  know  the  wavelet,  and  therefore  we  can  use  the  wavelet 
deconvolved  by  the  complex  cepstrum  as  the  input  signal  of  the  filter.  We  derive 
two  different  but  complementary  filters:  a  causal  Wiener  filter  which  assumes  a 
minimum  phase  wavelet  and  am  anti-causal  filter  which  assumes  a  maximum  phase 
wavelet.  This  section  considers  (a)  the  assumptions  about  the  signals  and  (b)  the 
derivation  of  the  Wiener  filter. 


4.1.  ASSUMPTIONS  ABOUT  THE  SIGNALS 

We  recall  that  the  received  signal  has  the  following  form: 

s(m)  —  r(rn)  *  u>(m)  +  n(m)  —  y(m)  +  n(m), 


where 


y(m)  =  r(m)  *  u>(m). 


w(m)  and  r(m)  are  respectively  the  wavelet  and  the  medium  response.  The  station- 
arity  of  the  signals  and  the  noise  is  assumed  up  to  the  second  order.  We  know  the 
second-order  statistics  £(?^m)2]  and  £[r(m)s(m)],  or  equivalent  statistics— as  we 
will  3ee  further  on.  Under  these  assumptions  we  estimate  the  medium  response  by 
using  an  estimator  which  is  „  linear  function  of  the  observation  s(m)  and  is  given 

by 

f(m)  =  s(m)  *  fi(m). 


The  Wiener  filter  characterized  by  fi(m)  is  defined  by  the  minimization  of  the  mean 
quadratic  error 

/[  +~ 

e2  ~  E  I  |  ^  h(k)s{m  -  k)  -  r(m) 

\  lk=-oo 

where  E  indicates  the  expected  value.  The  filter  h(m)  is  characterized  by  its  length 
P  and  the  located  interval  [-L,  P  -  L  • f  1]  :  h{m\F,L).  P  and  L  are  respectively 
called  the  order  and  the  lag  of  the  filter. 
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4.3.  DERIVATION  OF  THE  WIENER  FILTER 


s  4.2  1.  Derivation  of  the  zero-lag  Wiener  filter  for  minimum  phase  signals 

The  zero-lag  Wiener  filter  characterized  by  h(m,  P)  is  defined  by  the  minimization 
of  the  mean  quadratic  error 


e1  =  E 


h(k)s(m  -  k)  -  r(m ) 


The  filter  h(m,  P)  is  characterized  by  its  length  P  and  the  located  interval  [0,  P  —  1], 


Derivation  of  the  normal  equations  Let  us  define  the  prediction  error  ei  by 


P-i  P-i 

ei  -  ^  h(k)y(m  -  k)  +  h(k)n(m  -  k)  -  r(m),  m  £  [0,  P  -  l], 

Jc  =  0  k  =  0 

Then  we  have 

e7  =  E  (efet) . 

Let  H  be  the  matrix 

[h(0),h(l), .  ...h(P  -  1 )), 

Y  the  matrix 

[y(m),y(m-  l),...,y(m-  P+  1)], 

and  N  the  matrix 

(n(m),  n(m  -  1), . . , ,  n(m  -  P  -f  1)] 

The  mean  quadratic  error  can  be  rewritten  in  the  following  form  for  each  m: 

e2  =  E  ([ YHt  +  NHt  -  r(m)}T[YHT  +  N HT  -  r(m)]) . 

If  we  assume  that  the  sequences  y(m)  and  r(m)  are  uncorrelated  with  the  noise 
n(m)  and  if  one  expands  the  r.h.s.  of  this  equality  one  gets 

e2  =  HE\YtY]Ht  +  HE[NtN)Ht  -  2 E\r(m)Y)HT  +  £[r(m)2), 


Rvv,  Aqn  and  Rry  are  respectively  the  autocorrelation  matrix  of  the  sequence  y(m), 
the  autocorrelation  matrix  of  the  noise  n(m),  and  the  correlation  vector  of  the 
sequence  y(m)  with  the  scalar  r(m).  We  denote  Rvy(k),  Rn,,(k)  and  Rrv(k)  respec¬ 
tively  the  quantities  £[y(m  +  /)y(m+/  +  A)]t  J?[n(m  +  /)n(m  +  f-f  *)),  E(r(m)y(m  +  k)]. 
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If  we  assume  a  white  noise,  Rnn  is  a  diagonal  matrix  £„„(())./,  where  /  is  the  unit 
matrix  and  the  mean  quadratic  error  e1  is  given  by 

e*  =  HRvyHT  +  Rnn(0)HHT  -  2 RryHT  +  £[r(m)2]. 


Now  we  want  to  minimize  the  mean  quadratic  error  and,  classically,  e*  has  a  global 
minimum  if  the  two  following  conditions  are  fulfilled: 

(1)  Vat*  =  0, 

(2)  Vft(V//c2)  's  positive  definite  , 

where  V and  \7«(VHeJ)  respectively  the  gradient  and  the  hessian  of  the 
mean  quadratic  error  e2 . 

The  first  condition  leads  to 

2 HRVV  +  2Rnn(0)H  -  2Rrv  =  0 

and  the  second  condition  requires  the  matrix  Ruv  +  f?.nn(0)/  to  be  positive  definite. 
Thus  the  minimum  is  reached  when 

2h(Rvv-h  Rnn(0)I)=  Rrv.  (11) 


The  linear  system  (11)  is  called  the  set  of  normal  equations  for  the  Wiener  filter 
h(m,  P),  and  explicitly  the  set  of  normal  equations  is 


p-i 

v  h(k,  P)[Rw(m  -  k)  +  P„„(0)i5(m  -  1?)]  =  Prj,(m), 
h=o 


me|0,P-l).  (12) 


Remark  I  if  the  matrix  Rvv  +  iinn(0)/  is  positive  definite,  one  can  directly  find 
the  solution,  and  e2  can  be  expanded  into  a  quadratic  form  as  follows: 

e2  =  [h  -  Rry{Rvv  +  P„n(0)J)-1 1  [pvv  +  Pn„(0'/]  [ HT  -  (. R ,  +  Pnn(0)/)-l^v] 
-  Ryr(Ryy  +  Rnn(0)I)Rlr  +  ^Hm)*]. 

e2  vanishes  if  and  only  if  ff-Prv(flvv  +  f2n„(0)/)_1  vanishes  and  the  filter  coefficients 
are  given  by  the  exact  solution 

H  =  Rry(Ryy+  Rnn(0)I) 
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In  Eq.  (12)  we  cannot  access  the  correlation  matrices  Rn  and  Rry ,  and  so  we  replace 
them  by  well-known  second-order  statistics. 


If  we  take  the  i-transform  of  both  sides  of  Eq.  (12),  we  get  after  some  straightforward 
derivations 


p- 1 

V  h{k,P)z-k 

fc=0 


£  (Rn(m)+  Rnn{Q)6(m))z 


L  m  =  -  oo 


E 


Rry{m)z  m. 


According  to  the  Blackman  and  Tukey  definition  of  the  Power  Spectral  Density,  one 
gets 

p- 1 

£  P)*-k  (TW(Z)  +  Rnn( 0))  =  rr„(2),  (13) 

*=0 

where  rw(z)  is  the  power  spectral  density  of  the  sequence  y(n)  and  rry(r)  is  the 
cross-spectral  density  of  the  sequence  r(n)  with  the  sequence  y(n). 

If  the  sequence  r(n)  is  uncorrelated 

rry(z)  =  rrr(r)iv(z)*  =  jerr(o)fV(r)*,  (h) 

and 

rw(i)  =  rrr(2)ru,u,(r)  =  R„(o)rww(z),  (is) 

where  rrr(z)  and  rwur(2)  are  respectively  the  power  spectral  density  of  the  sequences 
r(n),  to(n)  and  Rrr  the  correlation  function  of  r(n). 

Then  Eq.  (13)  can  he  rewritten 


E 


*=o 


(16) 


Coming  back  in  the  time  domain,  Eq.  (16)  assumes  the  form 
£  h(k,P)  (^Rwv,(m  -  fc)  +  <5(m  -  k)j  =  u>(-m),  m  €  [0,P-  1).  (17) 


In  this  new  set  it  does  not  matter  if  we  do  not  know  the  autocorrelation  Rww  exactly 
as  we  can  estimate  it — through  u»(f)  being  the  wavelet  deconvolved  by  the  complex 
cepstium.  However  we  do  not  know  R„„(0)/ Rrr(0),  and  so  we  have  to  estimate  it. 
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Estimation  of  Rnn(0)/Rrr(0)  The  principle  of  this  estimation  is  the  eigenvalue 
decomposition  of  the  correlation  matrices  Rw  +  Rnn(0)J  [26]  and  Rwu/.  According 
to  classical  linear  algebra  these  matrices  can  be  decomposed  into  the  following  forms: 

RUV  +  Rnn(0)I  =  UyEyVf, 

in  which  Ey  is  the  diagonal  matrix  diag  (<r,  ,  <rlt . . .  ,  <rP),  with  <r*  the  eigenvalues  of 

£„„  +  £n„(0)/. 

Rww  —  UwT,wu£  , 

Eu,  is  the  diagonal  matrix  diag  (at ,  aj, . . . ,  op),  where  the  a<  are  the  eigenvalues 
of  Rww.  The  columns  of  Uw  are  the  orthonormalized  eigenvectors  associated  with 
these  eigenvalues.  In  our  case  the  eigenvalues  are  given  by 

+  Rnn{ 0),  *  €  (l,/>], 

with  A,  the  eigenvalues  of  the  correlation  matrix  Rvv.  At  this  point  let  us  assume 
that  the  rank  of  Rvv  is  smaller  than  P.  Because  of  the  equality  (15)  the  eigenvalues 
A,  and  are  linked  together  by  the  relation 

A;  =  Apr(0)on. 

Thus  the  rank  of  the  matrix  Rvv  is  equal  to  the  rank  of  the  matrix  RWV).  The  previous 
assumption  is  equivalent  to  assuming  than  the  order  of  the  wavelet  is  smaller  than 
P.  On  the  basis  of  this  assumption  the  estimation  procedure  is  as  follows: 

(1)  estimate  the  correlation  matrix  Rvv  +  Rnn( 0);  compute  the  eigenvalues  Oi 
and  the  eigenvectors  of  this  matrix; 

(2)  estimate  the  correlation  matrix  compute  the  eigenvalues  a,  and  the 

eigenvectors  of  this  matrix; 

(3)  estimate  the  rsmk  Q  of  the  correlation  matrix  Rvv\ 

(4)  estimate  Rnn( 0)  by  taking  the  average  of  the  P  -  Q  smallest  eigenvalues  <r< 
(the  eigenvalues  are  arranged  in  decreasing  order): 

*nn(0)  =  ( P~Q )  £ 

'  '  <=9+i 

(5)  estimate  the  eigenvalues  A <  from  A<  =  <r<  -  flno(0). 

(6)  estimate  fZrr(0)  by  means  of  rtrt.(0)  -  Q ~l  A^/oj. 

(7)  compute  £„„(0)/#rr(0). 

The  rank  Q  is  estimated  by  applying  the  AIC  Akaike  test  to  the  correlation  matrix 
RWw  +  ^nr»(0)/f?rr(0)/  (15).  One  recalls  that  this  test  consists  of  estimating  the 
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order  of  a  model  at  the  minimum  of  the  function 

^  (number  of  free  parameters) 

/(«)  =  -*(«)  +  1 - jjT* - - 

where  #(9)  is  the  Maximum  Likelihood  Function  of  the  order  q  ,  N  is  the  number 
of  observations  and  P  the  order  of  the  correlation  matrix.  In  our  case  the  function 
f(q)  is  (see  Appendix  C) 


f(q)  =  (P-q)  In 


P 

*  >=S+1 


P 


M*i) 

i=(+l 


+ 


q(P  +  1  -  1/2  9) 

N 


(18) 


Figure  29  represents  the  eigenvalues  of  the  correlation  matrix  Rww  for  different 
transmitted  pulses  of  16  time-samples  length  (in  these  simulations,  the  wavelet  was 
exactly  the  transmitted  pulse): 

•  CW  signal  windowed  by  a  rectangular  window, 

•  CW  signal  windowed  by  a  rectangular  window, 

•  CW  signal  windowed  by  a  half-cycle  sinusoidal  window, 

•  CW  signal  windowed  by  a  Hanning  window, 

•  CvV  signal  windowed  by  a  Hamming  window. 

In  Fig.  30  we  present  the  eigenvalues  of  fZu;u,(0)  for  the  same  pulses  but  of  64  time- 
samples  length. 

In  Figs.  31  and  32  we  present  the  Akaike  functions  f(q)  applied  to  the  matrix 
Ryy  +  i?n„(0)/  for  pulse  lengths  of  16  and  64  time-samples  respectively. 

In  Figs.  33  and  34  we  present  the  estimate  of  Rnn(0)/  fi,.r(0)  for  the  four  windows 
mentioned  above  for  respectively  16  and  64  time-samples  length.  These  results  shows 
that  the  Hanning  window  is  the  one  which  is  best  at  discriminating  the  eigenvalues 
corresponding  to  the  wavelet  and  the  eigenvalues  corresponding  to  the  noise. 


Solution  of  the  normal  equations  Taking  into  account  the  Toeplitz  form  of  Rww  + 
12»m(0)/f?rr(0),  the  normal  Eqs.  (17)  are  solved  by  the  Levinson  Algorithm  (16). 

Stability  of  the  filter  Since  the  wavelet  u/(n)  is  minimum  phae,  Eq.  (15)  ensures 
the  stability  of  the  Wiener  filter  defined  by  h(k,  P). 
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■  4.2.2.  Derivation  of  the  sero-lag  Wiener  filter  for  maximum  phase  signals 

The  derivation  is  similar  to  the  derivation  for  minimum  phase  signals.  This  time,  the 
filter  h(m,  P)  is  characterised  by  its  length  P  and  the  located  interval  [-P  +  1,0). 


Derivation  of  the  normal  equations  Here,  the  prediction  error  e,  is  defined  by 


-p+i  -P+i 

£  h(k)y(m  -  k)  +  ^  h(k)n(m  -  k)  -  r(m)t  m  €[-/*  + 1,0) 
*=o  s=o 


and  we  want  to  minimise  the  mean  quadratic  error 


e  =  £(elei). 


Using  the  same  derivations  that  for  minimum  phase  signals,  we  gvst  to  the  set  of 
normal  equations  which  define  the  Wiener  filter  h(m,  P) 

-p+i 

]T  h(k,P)[Rvv(m-k)  +  Rnn(0)6(m~k))  =  Rrv(m),  m  6  (-P  +  1,0).  (19) 
*=o 


If  we  assume  that  the  sequence  r(n)  is  uncorrelated  and  by  means  of  derivations 
similar  at  the  minimum  phase  case,  Eq.  (17)  becomes 


-p+i  . 

£  h(k,P)(R% 

fc=0  ' 


„{k  -  m)  +  -'-I'-’ 6{k  -  m)j  -  to(-m), 


*nn(0), 

^rr(O) 


me|-P  +  1,0) 


(20) 


The  estimation  of  Ann(0)/Arr(0)  is  identical  to  the  estimation  for  minimum  phase 
signals.  Since  RWVI  +  Rn„{0)/ Rrr(0)I  has  a  toeplits  form,  the  solutions  of  Eq.  (20) 
are  obtained  by  the  Generalised  Levinson  Algorithm.  Since  the  wavelet  w(n)  is 
maximum  phase,  Eq.  (20)  ensures  the  stability  of  the  Wiener  filter  defined  by  h(k,  P). 
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5.  Combination  of  homomorphic 
deconvolution  and  Wiener  filtering 

As  we  saw  in  Sect.  4,  the  Wiener  filter  is  well  defined  for  a  minimum  or  maximum 
phase  input  sequence,  but  it  is  rather  unstable  for  a  mixed  phase  sequence.  Since 
the  received  signal  and  the  wavelet  are  mixed  phase  in  real  life  (see  Appendix  B), 
and  idea  is  to  factorised  the  received  signal  and  the  wavelet  into  their  minimum  and 
maximum  phase  components.  Then,  in  order  to  improve  the  deconvolution  method, 
we  can  apply  a  sero-iag  causal  Wiener  filter  to  the  minimum  phase  component  and 
an  sero-lag  anticausal  Wiener  filter  to  the  maximum  phase  component.  This  idea 
has  already  been  used  by  Oppenheim  et  al.  (10),  using  a  linear  predictor  instead  of 
a  Wiener  filter.  (Note  that  we  do  not  present  any  results  here  on  the  combination 
of  homomorphic  deconvolution  and  Wiener  filtering.) 


6.1.  FACTORIZATION  OF  THE  MIXED  PHASE  SIGNALS 

Let  us  assume,  as  in  the  previous  chapters,  that  the  received  signal  p(t)  is  the 
convolution  of  the  wavelet  u»(t)  with  the  medium  response  r(t) 

V(0  =  w(t)  *  r(t), 


or  in  the  r-dom&in 

K(z)  =  W(z)R(z). 

If  we  assume  that  y(z)  is  a  rational  transfer  function,  W(z)  can  be  factorised  as 
follows: 

W(i)  =  Wwta  P(z)WmM  p(z), 

where  Wmi,,  p(z)  and  p(z)  are  respectively  the  minimum  and  maximum  phase 
components  of  W'(z).  In  the  same  way,  R(z)  can  be  factorised  as  follows: 

R(z)  =  Amin  p(z)R  m%x  p(*)i 

where  p(z)  and  p(z)  are  respectively  the  minimum  and  maximum  phase 
components  of  -ft(z).  Therefore,  V(z)  cam  be  rewritten  in  the  following  form: 

K(z)  =  [Wmin  p(z)R„.iB  p(z)][Wra„  p(z)Rm„  p(z)|, 

or 

y(z)^  vml.  p(z)ymup(z), 
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^mln  p(*)  —  ^mln  j>(*)ffcmln  p(*)> 

Vm«  p(z)  =  wm„  p(z)R  m*x  p(*)- 
In  the  cepstral  domain,  the  previous  equations  become 

y(0  =  y min  p(0  +  jmu  p(t)> 

=  [Wmln  p(0  "b  ^min  p(0)  "b  1*^  mu  p(0  +  ^m»x  p(01 


According  to  the  properties  of  the  complex  cepstrum  recalled  in  Sect.  3,  ym|„  p(t)  i» 
equal  to  zero  for  the  negative  frequencies  and  ymtz  p(t)  is  equal  to  sero  for  the  posi¬ 
tive  frequencies.  Then,  by  applying  the  complex  cepstrum,  we  are  able  to  factorised 
y(t)  and  ui(t)  into  their  minimum  and  maximum  phase  components. 


5.1.  PROCEDURE  TO  DECONVOLVE  THE  WAVELET  AND  THE  MEDIUM 
RESPONSE 

We  first  apply  the  complex  cepstrum  to  the  received  signal  y(t).  We  filter  the  com¬ 
plex  cepstrum  y(t)  by  means  of  two  rectangular  windows.  The  first  window  is  defined 
for  the  positive  frequencies  in  order  to  extract  the  cepstrum  ymi„  p(t).  The  second 
window  is  defined  for  the  negative  frequencies  in  order  to  extract  ymu  p.  Then,  we 
low-pass  filter  ymia  p(t)  to  separate  wmin  p(t)  said  fmln  p(t),  and  we  high-pass  filter 
p(0  to  separate  wmmx  p{t)  and  rmM  p(t).  Thus  we  get  both  the  minimiun  and 
maximum  phase  components  of  the  received  signal  and  the  deconvolved  wavelet. 
The  next  step  of  the  procedure  consists  of  simultaneously  applying  a  causal  Wiener 
filter  hm i„  (t)  to  ymin  p(t)  with  u>m|0  p(t)  in  input  and  an  anti-causal  Wiener  filter 
(0  to  ymu  p(t)  with  p(t)  in  input. 

The  medium  response  r(f)  is  estimated  by  the  inverse  filter  hmin(t)  *  hm**  (t)  as 
follow: 

^(t)  =  hmia  (t)  *  hmu  (0  *  y(t)' 

The  deconvolution  procedure  is  shown  schematically  in  Fig.  8b.  The  system  D.  is 
nothing  other  than  the  complex  cepstrum  as  defined  in  Sect.  3. 
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6.  Results 


S.l.  RESULTS  OBTAINED  BY  SIMULATIONS 

The  results  presented  in  this  section  relate  to  active  sonar  reverberation.  They  are 
summary-type  representative  results,  proving  the  feasibility  of  the  methods,  but  also 
pointing  out  their  limitations.  The  transmitted  signals  are  windowed  CW  pulses  and 
the  reverberated  signals  are  received  on  a  horisontal  towed  array. 


■  6.1.1.  Reverberation  in  active  sonar 

We  assume  that  reverberation  is  measured  in  deep  water  with  a  low-frequency  om¬ 
nidirectional  source  and  a  towed,  horisontal  array.  We  are  looking  at  reverberated 
signals  after  be&mforming.  The  simulations  try  to  be  an  accurate  copy  of  the  exper¬ 
iments  carried  out  for  backscattering  studies  in  active  sonar  (surface,  volume  and 
bottom-layer  backscattering). 

The  simulations  are  described  in  terms  of  two  models. 


■  6.1.1a.  First  model 

The  scenario  is  depicted  in  Fig.  35.  Remember  that  it  is  the  simulated  signals  after 
beamforming  that  are  simulated.  The  simulation  does  not  take  any  beamforming 
processing  into  account. 


Description  of  the  signals 


Transmitted  signal  The  transmitted  signal  is  a  Hanning-windowed  CW  pulse.  The 
pulse  length  is  taken  as  a  parameter  of  the  simulation.  The  sampled  CW  pulse  is 
modelled  in  time  es 

i*(n)  =  sin(2x/o(n  -  1)](1  -  cos(2*((n  -  1)/L))], 

where  L  is  the  pulse  length  expressed  in  time-samp’es  and  /o  is  the  normalised 
frequency  of  the  CW  signal.  The  CW  pulse  is  represented  in  Fig.  36. 


Medium  impulse  response  The  reverberation  model  has  three  paths:  the  direct 
path,  the  reflection  on  the  surface  and  the  reflection  on  the  bottom.  We  do  not 
model  the  transfer  function  of  the  surface,  nor  the  bottom  transfer  function.  The 
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travel  times  along  the  three  paths  are  three  parameters  of  the  simulation,  and  de¬ 
pend  on  the  source  depth  and  its  distance  from  the  array,  and  the  array  and  wa¬ 
ter  depths.  td,ts, Ti  are  respectively  the  arrival  times  for  the  direct  path,  the 
surface  reflected  path,  and  the  bottom-layer  reflected  path.  The  three  discrete¬ 
time  paths  are  modelled  as  three  Dirac  at  the  time-samples  np,r%s  and  nowhere 
rD  =  no  At,  t$  =  r»s  At,  tl  —  At;  At  is  the  time  sampling  interval)  and  is  given 
by 

r(n)  =  -r$(n  -  nD )  +  rl$(n  -  ns)  -  r3$(n  -  rn.) 


Additive  noise  The  noise  is  charact  rized  by  its  spectrum  and  the  signal-to-noise 
ratio,  and  is  defined  as  the  response  of  a  linear  Alter  to  an  input  white  gaussian  noise 
(random  normal  sequence).  Because  of  the  frequency  step  feature  of  the  complex 
cepstrum,  we  are  interested  in  the  SNR  at  each  frequency.  Hence,  three  SNRs  are 
defined.  One  of  these,  called  SNRT,  is  the  transmitted  signal-to-noise  ratio  and 
another,  called  SNRR,  is  the  received  signal-to-noise  ratio.  These  two  SNRs  are 
processed  in  the  full  frequency  band  as  follows: 


SNRTdB  =  10  log 


SNRRdB  =  10  log 


/7oB  !*(/)!’ d/\ 
UB|A'(/)lJd// 
(7oBlS(/)|2d/\ 


where  X(f),  S(f),  and  N{f)  are  respectively  the  spectrum  of  the  transmitted  signal 
x(t),  the  received  signal  s(t),  and  the  noise  n(t);  B  is  the  frequency  band.  We  recall 
that  the  normalised  Hanning  CW  pulse  bandwith  is  given  by  the  well-known  relation 


Sew  -  4ir /£. 


The  third  signal-to-noise  ratio,  called  SNRF,  is  defined  at  each  frequency  of  the 
filtered  bandwith  as  follows: 


SNRFdB  = 


Description  of  the  processing  An  observation  time  of  256  time-samples  has  been 
used  and  the  pulse  length  is  equal  to  64  time-samples.  The  normalized  frequency 
/o  of  the  CW'  signal  is  equal  to  0.25.  First  of  all  the  received  signal  ia  band-pass 
filtered  in  frequency  with  a  rectangular  window  defined  by  the  lowest  normalised 
frequency  /min  and  the  highest  normalized  frequency  /mM.  These  two  frequencies 
are  given  here  as  /min  =  0.222  and  /mM  -  0.277.  Then  we  apply  the  band-pass 
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mapping  and  the  complex  cepstrum  to  get  the  deconvolved  wavelet.  The  medium 
impulse  response  is  deconvolved  by  both  complex  cepstrum  and  Wiener  filtering. 


Results  an d  their  interpretation  For  a  given  pulse  length  and  multipath  configu- 
ration,  we  first  look  at  the  effect  of  the  noise  on  the  deconvolved  signal  accuracy. 
The  results  are  summarised  in  Table  1  (the  pulse  length  is  64  time-samples,  the 
delays  td,ts,ti  are  respectively  equal  to  80,  110  and  170  time- samples). 


Table  1 


Classification  of  the  results  with  the  number  of  the  corresponding  figure.  Active  sonar 
simulation:  Hanning  CW  pulse  with  3  multiples 


SNRR 

(dB) 

Received 

signal 

Power 

spectrum 

Deconvolved 

wavelet 

Deconvolved 

medium 

response 

(cepstrum) 

Deconvolved 

medium 

response 

(Wiener) 

13 

Fig.  38 

Fig.  39 

Fig.  40 

Fig.  41 

Fig.  42 

8 

Fig.  43 

Fig.  44 

Fig.  45 

failed 

Fig.  46 

3 

Fig.  47 

Fig.  48 

Fig.  49 

failed 

Fig.  50 

-  1 

Fig.  51 

Fig.  52 

Fig.  53 

failed 

failed 

-  6 

Fig.  54 

Fig.  55 

Fig.  56 

failed 

failed 

-  11 

Fig.  57 

Fig.  58 

Fig.  59 

failed 

failed 

We  conclude  that  the  wavelet  is  rather  well  deconvolved  up  to  a  SNRR  of  -6  dB 
and  seems  relatively  insensitive  to  additive  noise.  The  wavelet  can  be  rescaled  by 
correlation  with  the  transmitted  pulse.  The  correlatir  a  function  of  the  deconvolved 
wavelet  with  the  CW  pulse  for  a  received  signal- to- noise  ratio  of  16  dB  is  depicted 
in  Fig.  37.  On  the  other  hand,  the  medium  impulse  response  suffers  more  from  the 
additive  noise,  as  we  can  see  in  Fig.  41.  The  complex  cepstrum  cannot  deconvolve 
the  medium  response  at  lower  signed- to-noise  ratio.  These  results  agree  with  the 
mathematical  derivation  in  Appendix  B,  where  it  is  shown  that  the  medium  impulse 
response  deconvolved  by  the  complex  cepstrum  is  more  affected  by  additive  noise 
than  the  wavelet  is.  The  Wiener  filter,  with  the  original  pulse  as  input,  can  accept¬ 
ably  separate,  the  three  multipaths  shown  in  Fig.  42.  It  does  so  successfully  up  to  a 
SNR  of  3  dB  (the  corresponding  SNRF  values  are  given  in  Table  2).  In  Figs.  42,  46 
and  50,  we  see  that  the  Wiener  filter  resolution  is  not  ‘optimal’,  as  a  consequence  of 
the  fact  that  we  add  a  ‘white  noise’  parameter  to  the  zero-lag  element  of  the  auto¬ 
correlation  matrix  in  order  to  stabilize  the  computation  of  the  inverse  [24j.  Here,  the 
white  noise  parameter  is  equal  to  0.005.  The  ill-conditioned  problem  arises  because 
the  order  of  the  received  signal  is  smaller  than  the  order  of  the  linear  system,  as 


-31- 


Saclantcbn  SM-aos 


shown  by  the  Akaike  test  in  the  Sect.  4.  Another  reason  is  the  non-minimum  phase 
characteristic  of  the  transmitted  pulse. 


Table  2 

Signal-to-noise  ratio  of  the  received  signal  at  each 


frequency  of  the  bandwidth: 
with  3  multiples 

Hanning  CW  pulse 

Normalised  frequency 

SNRF 

(dB) 

0.222 

-  90 

0.220 

6 

0.230 

8 

0.234 

4 

0.238 

4 

0.242 

12 

0.246 

17 

0.250 

18 

0.254 

17 

0.258 

24 

0.281 

17 

0.285 

14 

0  289 

12 

0.273 

-  9 

0.277 

-  4 

SNRT  =  10  dB  SNRR  =  3  dB 

■  6.1.1b.  Second  model 


Description  of  the  signals 


I Yansmitted  signs’  The  same  transmitted  pulse  described  in  the  previous  section 
is  used. 

Medium  impulse  response  Five  equi-spaced  multiples  defined  by  the  time  delays 
’I.  v«,  and  n,  with  values,  respectively,  of  80,  110,  140,  170  and  200  time- 

samples. 
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Additive  noise  The  same  characteristics  as  that  in  the  previous  section. 


Description  of  the  processing  The  same  processing  as  that  in  the  previous  section. 


Results  and  their  interpretation  The  results  are  summarized  in  Table  3. 


Table  3 


Sign&l-to-noise  ratio  of  the  received  signal  at  each  frequency  of  the  bandwidth:  Hanning 
CW  pulse  with  3  multiples 


SNRR 

(dB) 

Received 

signal 

Power 
sf  ?ctrum 

Deconvolved 

wavelet 

Deconvolved 

medium 

response 

(cepstrum) 

Deconvolved 

medium 

response 

(Wiener) 

16 

Fig.  61 

Fig.  62 

Fig.  63 

Fig.  64 

Fig.  65 

11 

Fig.  66 

Fig.  67 

Fig.  68 

Fig.  69 

Fig.  70 

6 

Fig.  71 

Fig.  72 

Fig.  73 

Fig.  74 

Fig.  75 

1 

Fig.  76 

Fig.  77 

Fig.  78 

Fig.  79 

Fig.  80 

-  4 

Fig.  81 

Fig.  82 

Fig.  83 

failed 

failed 

Table  4 

Location  of  the  poles  of  jYo(z)  given  in  polar  coordinates 


Pole 

(no.) 

Radius 

Angle 

(dg) 

Pole 

(no.) 

Radius 

Angle 

(dg) 

1 

0.1 

0 

11 

0.7 

2 

0.2 

10 

12 

0.75 

110 

3 

0.3 

20 

13 

0.8 

120 

4 

0.3 

30 

14 

0.8 

130 

5 

0.4 

40 

15 

0.8 

140 

6 

0.45 

50 

16 

0.85 

150 

7 

0.5 

60 

17 

C.85 

160 

8 

0.55 

70 

18 

0.9 

170 

9 

0.6 

80 

19 

0.95 

180 

10 

0.65 

90 

20 

0.5 

180 
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We  conclude  that  the  wavelet  is  well  deconvolved  up  to  a  SNRR  of  0  dB.  The 
deconvolved  wavelet  can  be  rescaled  by  correlation  with  the  transmitted  pulse  (as 
in  the  previous  model)  Figure  60  depicts  this  correlation  function  for  a  SNRR  of 
15  dB.  The  deconvolutions  of  the  medium  impulse  response  by  the  complex  cepstrum 
and  Wiener  filtering  are  equivalent,  and  successful  up  to  0  dB.  The  time  delay  in  the 
medium  response  deconvolved  by  the  complex  cepstrum  is  due  to  the  linear  phase 
not  being  recovered  properly,  but  the  relative  positions  of  the  multiples  are  correct. 


6.2.  RESULTS  OBTAINED  WITH  EXPERIMENTAL  DATA 

■  6.2.1.  Reverberation  in  active  sonar 


Experiment  configuration  The  purpose  of  the  experiment  was  to  measure  rever¬ 
beration  with  an  activated  towed  array  at  low  frequency.  The  geometrical  configu¬ 
ration  is  depicted  in  Fig.  84  (23).  The  towed  array  has  32  hydrophones  spaced  at 
one  hall- wavelength  (1.96  m  for  the  measurements  processed  here).  The  array  depth 
was  around  100  m  and  was  separated  from  the  towship  by  900  m.  The  water  depth 
was  around  3500  in. 

Signal  characteristics  The  transmitted  signal  was  a  Hanning-windowed  CW  pulse. 
Its  duration  was  2  s  at  a  frequency  of  370  kHz.  The  signal  received  on  the  array 
is  beamformed  and  band-pass  filtered  in  frequency.  The  sampling  frequency  at  the 
beamformer  output  was  70  Hz,  and  the  observation  time  was  3.65  s  (256  time- 
samples).  The  resolution  in  time  provided  by  this  transmitted  signal  is 

_ 1 

bandwidth  ’ 


or  in  this  case 


tr  =  1.438  s. 


The  blocks  of  recorded  data  are  characterized  by  the  number  of  the  beam,  the 
number  of  the  ping,  and  the  range. 

Description  of  the  processing  The  processing  was  the  same  as  that  for  the  simu¬ 
lated  data  of  Subsubsect.  6.1.1a. 

Results  Figure  85  depicts  the  transmitted  pulse  in  time.  The  received  signal 
corresponding  to  beam  6,  ping  9  and  range  8  is  represented  in  Fig.  86.  The  wavelet 
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deconvolved  by  means  of  the  complex  cepstrum  is  represented  in  Fig.  87.  When  we 
compared  the  deconvolved  medium  response  for  two  different  inputs  of  the  Wiener 
filter  (the  deconvolved  wavelet  and  the  transmitted  pulse).  The  results  depicted 
Figs.  88  and  89  were  obtained.  The  results  are  similar  in  both  figures,  except  that 
we  get  a  better  resolution  with  the  transmitted  pulse  (Fig.  89  The  deconvolved 
wavelet  looks  like  the  transmitted  pulse,  which  is  a  promising  i  suit.  However,  the 
pulse  length  reduces  the  credibility  of  the  results  significantly,  j  •  bviously  not  the 
type  of  signed  one  should  use  to  study  reverberation.  (In  order  to  measure  surface, 
volume  and  bottom  backscattering,  the  transmitted  pulse  must  have  a  significant 
bandwidth,  and  instead  of  using  CW  pulses  of  2  s  it  would  be  more  sensible  to 
transmit  pulses  of  0.1  s,  for  example.) 


■  6.2.2.  Explosive  data 


Experiment  configuration  The  data  come  from  an  acoustic  propagation  experi¬ 
ment  made  by  the  Centre’s  Environmental  Acoustics  Group  in  the  Tyrrhenian  sea. 
The  aim  of  the  experiment  was  to  estimate  the  transfer  function  of  the  ocean  over 
a  broad  acoustic  frequency  range.  The  broadband  signal  arising  from  an  explosive 
source  was  recorded  (a)  at  a  range  of  4.5  km  with  a  vertical  array  of  32  hydrophones 
spaced  at  2  m,  and  (b)  close  to  the  source  with  a  portable  array  of  4  hydrophones. 
The  experiment  configuration,  with  its  various  geometrical  parameters,  is  presented 
in  Fig.  90.  Before  any  kind  of  processing  one  can  expect  at  least  four  arrivals:  one 
direct,  one  by  reflection  and  one  by  refraction  at  the  sea  surface,  and  later  one  by 
reflection  at  the  seabed.  Since  it  is  a  deep  water  environment,  we  do  not  consider 
the  bottom  reflection. 


Signal  characteristics  The  explosive  is  a  broad  band  source.  The  power  spectrum 
of  the  signal  received  on  hydrophone  17  of  the  vertical  array  is  depicted  in  Fig.  93. 
The  sampling  frequency  was  6  kHz.  Therefore,  according  to  Fig.  93  the  received  sig¬ 
nal  band*ith  was  almost  2.5  kHz.  The  time  series  at  the  output  of  each  hydrophone 
is  represented  on  Fig.  91. 


Description  of  the  processing  The  observation  time  was  170  ms  (or  1024  time- 
samples).  We  processed  the  full  frequency  band  (no  band-pass  mapping).  The 
wavelet  was  deconvolved  by  the  complex  cepstrum,  and  the  paths  reflected  and 
paths  refracted  at  the  surface  were  resolved  by  Wiener  filtering. 


Results  The  received  time  series  for  hydrophones  4  and  17  are  depicted  on  Figs.  92 
and  94.  The  deconvolved  wavelets  are  presented  in  Figs.  95  and  96  The  minimum 
phase  property  of  the  wavelet  can  be  studied  in  Fig.  97,  which  shows  the  partial 
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energy  of  the  two  deconvolved  wavelets.  The  partial  energy  is  defined  as  follows  [22] 


£p(m)  =  Mn)|*. 


The  two  deconvolved  wavelets  carry  the  same  quantity  of  energy.  The  wavelet  for 
hydrophone  4  has  more  energy  at  the  beginning;  the  wavelet  for  hydrophone  17  has 
more  energy  at  the  end.  Their  power  spectra  are  identical,  as  we  can  see  in  Figs.  98 
and  99.  It  seems  that  some  poles  or  zeros  of  the  transfer  function  of  the  wavele  has 
been  transferred  outside  of  the  unit  circle.  We  ct\n  see  that  the  wavelets  are  not 
minimum  phase,  most  definitely  for  the  wavelet  corresponding  to  hydrophone  17.  If 
we  compare  the  results  to  a  theoretical  wavelet,  it  seems  the  original  shot  has  been 
perturbed  by  the  propagation  medium  and  perhaps  also  the  layer  conditions  (the 
sea-surface  was  flat  during  the  experiment  and  introduced  only  a  time-delay).  If  we 
use  the  deconvolved  wavelet  in  order  to  resolve  the  reflected  and  refracted  paths,  we 
have  no  success.  Therefore,  we  use  the  first  arrival  as  the  input  of  the  Wiener  filter. 
The  deconvolved  reflected  and  refracted  arrivals  corresponding  to  hydrophones  4 
and  17  are  presented  respectively  in  Figs.  100  and  101.  The  resultc  for  the  entire 
vertical  array  are  presented  in  Fij  .  102,  with  the  direct  arrival  taken  as  the  time 
origin. 

These  results  confirm  the  hypothesis  of  three  main  arrivals,  one  direct,  one  refracted 
and  one  reflected.  The  assumed  propagation  model  is  depicted  in  Fig.  102b,  which 
also  shows  the  mean  sound-velocity  profile  estimated  from  the  measurements. 


-36- 


SaCLANTCEN  SM-203 


7.  Conclusions 


The  study  carried  out  in  this  report  has  pointed  out  the  importance  of  the  phase 
information  in  the  understanding  of  propagation  and  reverberation  mechanisms. 
Since  the  phase  behaviour  is  rather  complicated  due,  in  part  to  its  randomness,  we 
need  some  accurate  signal  processing  methous  to  perform  the  analysis.  By  modeling 
the  propagation  medium,  the  bottom-layer,  and  the  surface  as  linear  filters  one  may 
apply  techniques  such  as  deconvolution  and  identification  methods. 


The  complex  cepstrum  used  to  deconvolve  the  wavelet  does  not  postulate  a  minimum 
or  maximum  phase  characteristic  for  the  signals,  and  therefore  is  very  useful  in 
propagation  and  reverberation  application,  for  which  the  signals  are  mostly  mixed 
phase.  Although  the  complex  cepstrum  requires  a  relatively  high  signal-to-noise 
ratio,  because  of  the  phase  unwrapping,  its  application  to  seismic  and  active  sonar 
reverberation  is  meaningful.  The  results  obtained  and  presented  in  Sect.  6  with 
simulated  data  are  quite  satisfactory  up  to  a  signal-to-noise  ratio  of  -5  dB  for  the 
wavelet,  but  only  up  to  5  dB  for  the  medium  response.  These  results  confirm  the 
derivation  made  in  Appendix  B,  proving  that  the  complex  cepstrum  does  not  succeed 
in  deconvolving  the  medium  response  as  well  as  it  does  for  the  wavelet,  due  to  the 
presence  of  additive  noise.  We  have  seen  that  the  deconvolved  wavelet  carries  a 
lot  of  information  on  the  medium  and  a  further  focus  would  be  to  fit  a  parametric 
model  and  control  its  behaviour  with  the  propagation  conditions.  This  can  be  done 
by  using  autoregressive  (AR)  or  autoregressive-moving  average  (ARMA)  modelling 
of  the  wavelet. 


The  results  shown  in  Sect.  6,  on  explosive  measured  data,  emphazise  the  impor¬ 
tance  of  the  phase.  We  have  seen  that  the  deconvolved  wavelet  at  two  separate 
hydrophones  of  the  vertical  array  can  have  the  same  power  spectrum  but  not  the 
same  energy  distribution  (shown  by  the  partial  energy  curves)  The  transfer  function 
treated  in  Appendix  A  is  an  example  of  how  an  all-pass  filter  can  modify  the  energy 
distribution  inside  a  signal,  transforming  it  from  a  minimum  phase  signal  to  a  mixed 
phase  signed.  This  can  arise  when  the  transmitted  pulse  goes  through  a  layer  which 
has  an  all-pass-filter  transfer  function  (Appendix  B).  In  this  environmental  technique 
like  the  power  cepstrum  method  is  not  powerful  enough  One  deficiency  underlined 
in  this  report,  concerning  the  complex  cepstrum  technique,  is  the  restoration  of  the 
linear  phase  once  the  different  components  have  been  deconvolved.  In  the  described 
simulations  in  Sect.  6,  tliis  equipment  has  been  solved  by  ‘omputing  the  crcsscorre- 
lation  between  the  deconvolved  wavelet  and  the  transmitted  pulse.  However,  some 
more  attention  has  to  be  put  on  this  particular  but  significant  part  of  the  complex 
cepstrum. 
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The  combination  of  homomorphic  deconvolution  and  Wiener  filtering  is  well  adapted 
to  reverberation  studies  in  that  it  capitalises  on  the  individusd  advantages  of  both 
the  techniques.  The  homomorphic  deconvolution  handles  the  mixed  phase  char¬ 
acteristics  of  the  wavelet  while  the  WieneT  filter  provides  high  resolution  of  the 
medium  response,  as,  for  example  of- arrivals  estimation.  Some  results  using 
this  promising  mixed  technique  will  be  the  subject  of  a  subsequent  report. 
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Appendix  A 

Minimum-phase  signals  and  their  properties 

The  notion  of  minimum-phase  signals  is  of  considerable  importance  in  signal  process¬ 
ing  [9].  In  this  section  it  will  be  shown  that  the  fourier  tranform  of  a  minimum-phase 
signal  can  be  recovered  from  its  magnitude  or  its  phase.  Because  most  of  the  digital 
Alters  are  defined  in  term  of  magnitude,  it  is  important  to  know  the  phase  in  order 
to  design  stable  filters.  The  minimum-phase  condition  gives  some  nice  properties  to 
the  complex  cepstrum  (8]  and  allow  us  to  design  inverse  filters  [16]. 

Before  giving  the  definition  of  a  minimum  phase  signal,  let  us  recall  the  definition 
of  a  causal  signal  in  order  to  make  an  analogy  between  the  complex  cepstrum  of  a 
minimum-phase  sequence  and  a  causal  sign&i. 


A.l.  DEFINITION  OF  A  CAUSAL  SIGNAL-PROPERTY  OF  ITS  FOURIER 
TRANSFORM 

The  values  of  a  causal  signal  *(t)  are  null  for  the  negative  values  of  t.  We  can  always 
write  *(f)  as  a  sum  of  an  odd  function  *odd(t)  and  even  function  : 

*(t)  =  *odd(<)  +  (0. 


where 


*.».„(*)  =  £(*(0  +  *(-<)),  (A.l) 

*.«(*)=  H*(«)  -*(-«)>•  (A. 2) 


For  the  positive  values  of  t  we  have 

*«*«n{0  =  *odd(0  = 


and  for  the  negative  values 


*odd(0  —  *tT«n(f)  —  —  t)i 


and  we  can  rewrite 

*odd(^)  =  (A. 3) 

*m«(0  =  sig(t)Xodd(0’  (A.  4) 
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By  taking  th«  Courier  transform  of  both  of  the  tides  of  the  equality  (Eq.  A.l)  and 
calling  respectively  A.T„,(/)  and  X{f)  the  Courier  transform  of  x,*«B(0  And  *(1)  , 
we  have 


*„..(/)  =  r~  dt  =  f  r°°  z(t)e-»"  .U  if  [  +  a°  ri-fy-^dt 

«/  —  OO  w  *  OP  J  ~OL> 

=  f  r  zwe-’*1'*  dt  -  \  r  x(t)cj,'‘d< 

«/  —  OO  w  —  OO 

=  f[A(/)-A'(/))  =  flmT(/). 


By  taking  the  Courier  transform  of  both  the  sided  of  Eq.  (A. 2)  and  calling  Aodd(/) 
the  Courier  transform  of  Co4d(t)<  we  get 

Aodd(/)=  f  /+0O*(t)e->*-*dt-  f  [+~  z(-t)e*Ut  dt 

J -OO  -OO 

=  f|A(/)  +  A(/)]  =  Re*(/). 


If  we  take  now  the  Courier  transform  of  (A. 3)  we  have 

Aodd(/)  =  Courier  transform(sig(t))  *  A.„0(/), 

=  —  vp(i)*tImA(/), 

where  vp(lj f)  is  the  Cauchy  principal  value  of  /**(l//)d/. 

Then 

v  1  /+0°  »1mA(v) 

Aodd (/)  =  —  vp  I  — - '-dv 

•*  /-»  /-*' 

or 

ReA(/)  =  -vp  /  1  dv  =  hilbert  transform(bnx(/)j. 

*  J -CO  I  ~~  v 

II  we  take  the  Courier  transform  of  (A. 4),  we  have 


A„.„(/)  =  ~vp(I)*ReA(/), 

KsX(u) 


.1  f+°°1i3X 

'«vpL  f- 


diJ 


=  inverse  hilbert  traasform(Re(A(/)j. 
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By  performing  the  inverse  process  we  show  that  if  Im  X{f)  and  Re X(f)  are  related 
through  a  hilbert  transform,  the  signal  z(t)  is  causal. 

Therefore  a  causal  signal  is  characterised  by  the  fact  that  the  real  part  and  the 
imaginary  part  of  its  fourier  transform  are  related  through  a  hilbert  transform. 

The  z-transform  of  a  causa]  sequence  z(n)  converges  in  the  domain 

(R,+ao) 

but  this  necessary  condition  is  rot  sufficient  to  steady  the  causality  of  z(n).  A 
necessary  and  sufficient  condition  is  given  by  the  following  theorem. 

Theorem:  A(z)  is  the  z-transform  of  a  causal  sequence  *(n)  if  and  only  if  2f(z)  is 
bounded  when  z  reaches  the  infinity. 


A  J  DEFINITION  OF  A  MINIMUM-PHASE  SIGNAL 


■  .4.2  1.  Definition 
Let  A(z)  be  defined  by 

X(z)  -  logjf(z)  =  log|AT(z)|  -f  *argA(z), 

and  let  i(n)  be  the  inverse  z-transform  of  A(r);  i(n)  is  by  definition  the  complex 
cepstrum  of  z(n)  [see  Sect  3]. 

The  minimum-phase  condition  is  that  the  complex  i  epstrum  i(n)  is  causal  or,  ac¬ 
cording  to  the  previous  section,  that 

lcg|A(z)|  =  hilbert  tranzfcrm(arg  JY(z)|.  (A  5) 

For  a  minimum-phase  signal  z{£),  the  phase  o.  X (/)  is  uniquely  defined  from  the 
magnitude  log  |A'(/)|.  Another  condition  is  that  there  is  causal,  stable  inverse  system 
X ~*(z)  such  that 

x-l(z)jr(z)  =  i 


There  is  a  consequential  a  property  of  the  minimum-phase  sequences  x(n)  the  poles 
and  the  seros  of  the  z  transform  X(z)  are  inride  the  un:t  circle. 


■  A. 2.2.  Justification  of  the  terminology  ‘minimum-phase  signal'  [22] 

Let  X(z)  be  the  z-transform  of  any  signal  *(t).  A(z)  can  be  written  as  a  product 
of  two  functions 


X(z)  =  Ao{z)(?*  p.(z), 


(A.3) 
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where  Xo(z)  is  the  z-transform  of  a  minimum-phase  signal  and  G»  p  (z)  the 
transfer  function  of  an  ail-pass  filter.  The  role  of  G,.p.(z)  is  to  transfer  M  zeros  of 
Xo(z)  outside  the  unit  circle  without  modifying  the  magnitude  response  |A0(z)|  and 
X(z)  can  be  rewritten  in  the  following  form 

M 

X(z)  =  G0(z)  J](z-  zm), 

msl 


where  zm  -vre  the  M  zeros  of  X{z)  outside  the  unit  circle.  To  render  the  terminology 
‘mini.Tium- phase’  explicit,  let  G,.p  (z)  be  an  all  all-pass  transfer  function  of  order 
one  w  th  only  one  real  zero.  G,  p  (z)  is  given  by 

G,  p  (z)  =  T~7~~  with|c|<l. 

14 -  cz 

G,  p  (z)  has  a  zero  (Z)  at  z(  =  -c  and  a  pole  (P)  at  zi  =  The  zero-pole  diagram 
of  G,  p  (z)  is  presented  in  Fig.  103. 

The  phaie-lag  angle  of  G.  p  (z)  is  given  by 

*(*)  =  -($;?(*)  -  */>(*))  =  *p(*)  - 

where  d'z(z)  and  #p(z)  are  respectively  the  angles  of  the  vectors  PA?  and  ZA} 
with  the  axis  0*.  For  the  normalized  frequency  /  in  the  range  [0,  the  phase-lag 
is  always  positive.  Let  $*(z)  and  $*0(z)  the  phases  of  X(z)  and  A<>(z),  and  we 
have 

**(*)  =  **„(*)-  *(*) 
or 

-$jc(z)  =  -#*0(z)  +  *(z),  (A. 7) 

and  therefore  the  phase-lag  angle  of  the  function  X(z)  is  always  greater  than  the 
phase-lag  angle  of  the  function  Ao(z)  The  all-pass  transfer  function  G».p.(z)  can 
be  decomposed  into  a  product  of  M  all-pass  transfer  functions  of  order  one,  and 
consequently  following  (A. 7)  the  phase-lag  angle  of  any  function  ^(z)  is  always 
greater  than  the  phase-lag  angle  of  the  function  ..Yo(z),  which  is  why  zo(0  is  called 
a  minimum  phase-lag  signal  or  by  abreviation  minimum-phase  signal. 

Remark  Given  two  functions  Jfi(z)  and  Aj(z),  one  has 

|A-,{z)|  =  |A,(z)|  and  #*,(z)  *  ♦*,(*) 
if 

*i(z)  =  Xo(*)G..Pl(*) 

A,(z)  =  Jfo(z)G,.p.,(z). 
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Hence  for  a  given  magnitude  |Jf(z)|  =  |Xi(z)|  =  |Ai(z)|,  the  phase- lag  is  not  defined 
uniquely.  It  is  uniquely  defined  if  *i(t)  and  x} (t)  are  minimum  phase. 


Now,  let  us  consider  examples  of  minimum-phase  and  mixed-phase  signals.  Let 
Jt0(z)  be  the  z-transform  of  a  minimum-phase  signal  x0 (t).  Aro(z)  is  represented  by 
an  all-pole  model  of  order  20  as  follows 

*°(2) = n  (7~Ki 1 


where  represents  the  poles  of  the  z-transform  located  inside  the  unit  circle.  Their 
positions  are  shown  in  Fig  104  and  the  exact  values  given  in  Table  4.  The  signal 
Zo(t)  is  represented  in  Fig.  107.  Let  xi(i)  be  the  signal  obtained  by  applying  an 
all-pass  filter  to  xo(t).  The  all-pass  filter  is  defined  by  its  transfer  function  in  the 
z-domain  by 


<w.)=n 

i=i 


(i  -PC*)' 


The  all-pass  filter  moves  the  seven  first  poles  of  Ao(z)  outside  the  unit  circle.  The 
poles  of  A'x(z)  are  represented  in  Fig.  105;  the  signal  Xx  (t)  is  represented  in  Fig.  108. 
Let  xj (t)  the  signal  obtained  by  applying  to  x0(f)  the  all-pass  filter  defined  by 


(0 


~  n 

i=i 


(‘-Pi) 

(i  -PC‘Y 


This  all-pass  filter  moves  the  ten  first  poles  of  Ho(z)  outside  of  the  unit  circle.  The 
poles  of  Xt(z)  are  represented  in  Fig.  106;  and  xj(t)  in  Fig.  109.  The  partial  energies 
of  the  three  signals  are  compared  in  Fig.  110.  We  recall  that  the  partial  energy  of  a 
signal  x(t)  is  defined  by 

m 

£(m )  =  53  *(«)*  • 

i  =  l 

The  minimum-phase  signal  is  the  one  which  has  the  energy  concentrated  at  the 
beginning 

Thus  the  all-pass  filter  can  partially  model  some  bottom-layers  and  incidentally  show 
their  influence  on  a  propagating  wavelet. 
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A  S.  RELATION  BETWEEN  THE  POWER  CEPSTRUM  AND  THE  COMPLEX 
CEPSTRUM 

We  recall  that  the  power  cepstrum  proposed  by  Bogert,  Healy  and  Tukey,  is  defined 

by 

z(t)  =  (fourier  t^an1fo^m(log|A'(/)|,)], 

where  X(f)  is  the  fourier  transform  of  r(t).  Because  X(f)  =  log  X(f)  =  log  |Ar(/)|  + 
» arg  X(f)  ~  XH{f)  +iX[(f)  we  have 

log  !AT(/)I*  =  2XH{f), 

because  Xn(f)  is  even  function  of  the  frequency  / 

—  oo 

rC°2XR(f)e»*'tdf\i. 

-oo 


The  integral  XR(f)ejiw,t  df  is  the  even  part  of  the  complex  cepstrum  x(t), 
denoted  by  i,v.n(t)  Consequently  we  have  the  relation 


One  can  always  decompose  the  complex  cepstrum  x(t)  to  the  sum  of  its  even  part 
and  its  odd  part 

i(t)  =  *,„„(<)  +  rodd(t), 

where  *„.„(<)  =  $(*(*)  +  *(“0)  “»d  *«dd(0  =  i(*(<)  -  *(-<))• 

If  the  complex  cepstrum  vanishes  for  the  negative  values  of  t, 

*.».n  (<)  = 


and  then 

!*(<)!’  =  r(t).  (A. 8) 

We  will  see  that  the  complex  cepstrum  vanishes  for  the  negative  values  of  t  in  the 
case  of  minimum  signals 

We  can  conclude  that  for  a  minimum  phase  input  sequence,  the  power  cepstrum  and 
the  complex  cepstrum  are  equivalent. 
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Appendix  B 

An  example  derivation  of  the  complex  cepstrum 

Jn  this  section  we  consider  a  two- path  propagation  model  with  reflection  at  the 
bottom  and  the  surface.  The  backscattering  mechanism  at  each  of  the  boundaries 
is  modelled  by  a  linear  filter.  The  objective  is  to  show  how  the  modelling  of  the 
propagation  medium  affects  the  complex  cepstrum  technique  (with  regard  to  the 
fundamental  notion  of  minimum  phase  signals)  and  to  point  out  the  limitation  of 
this  technique  for  deconvolution  of  the  medium  response.  The  propagation  model  is 
simple  and  highlights  the  parameters  which  have  a  determining  effect  on  the  method. 


B.l .  DERIVATION  OF  THE  COMPLEX  CEPSTRUM  FOR  THE  TWO-PATHS 
DISTORTION  CHANNEL  PLUS  NOISE 


a  B.1.1.  Expression  of  the  complex  cepstrum 

The  general  scenario  is  described  in  Fig.  111.  Here,  we  consider  the  case  where 
the  direct  path  is  not  taken  into  account.  In  practice,  this  means  that  the  observa¬ 
tion  time  starts  with  the  direct  arrival.  Under  this  assumption,  the  received  signal 
assumes  the  form 


•s(t)  =  h,(i  -  ri)  *  x(t)  +  h.(t  -  r.)  *  x(t)  +  n(<), 


where  we  recall  that  h((<)  and  h,(t)  are  respectively  the  impulse  response  associated 
with  the  bottom  layer  and  the  surface  and  the  additive  noise  n(t);  and  r,  are  the 
reflected  time  arrivals  with  respect  to  the  direct  arrival.  In  the  frequency  domain 
the  equality  becomes 


S(w)  =  +  H,(u)X{u>)e-i“T'  +  N{u>). 

Equation  (B.l)  can  be  factorized  into  the  following  form: 


(B.l) 


S(o>)  =  X(w)Ht(w)e~ivT‘ 


t  _L  ^ ,  N(w) 

+  H,(u)  +  Ht(u)X(w)  J 


By  normalizing  S(l>)  (removing  the  linear  phase  e  ■'wr)  and  taking  the  complex 
logarithm  of  both  jides  of  the  initial  equality),  we  have 


S(  u>)  =  X(u>)+  £,(*)  + log 


lI+ 


^  -jv(  r,-ri)  , 


1 

Ht(u,)X(u,)  i  ’ 
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where 


5(w)  -  logS(w),  =  log  A’(w),  =  log 


Under  the  assumption  that 

Hlju)  -juj-r.-Tt)  ,  N(u/)  JWTi  ^  , 

H,(u>)  Ki(u>)X(u>)  '  ' 

the  logarithm  can  be  expanded  into  its  Taylor  series  and  5(w)  takes  ti  e  form 


SH  =  X(u-)  +  Ht(u)  +  +  — yjffi  -< 


_  i  r  //«(w) 

*  Ht(u)) 


-ju2(t, -r,) 


H,{u>)  J 


’U( 


(^)A’(ur). 


//,(u;)JX(w) 


ftemarlc  1  The  medium  response  introduces  a  linear  phase  term  e  ‘■’u'Tiand  thete- 
fore  the  mean  phase  derivative  of  a(t)  is  not  eq  lal  to  0.  After  deconvolution  tnis 
Unear  phase  must  be  restored  to  the  deconvolves  medium  response. 

The  complex  cepstrurr  is  obtained  by  taking  the  inverse  fourier  transform  of  S(u>) 
and  is  given  by 


i(f)  =  i{t)  +  ht(t)  T  /»*[<  -  (r,  -  r,)]  +  hy(t  +  tj)  -  ljhi(t)  *  h2[t  -  2 (r,  -  r,)] 


-  *  *3(<  +  2 r,)  -  hi(t  - 


(B-2) 


where  fij(t),  fi3(r)  and  /»<(<)  are  respectively  the  inverse  fourier  transforms  of 


H.{u>)  N(u)  H,(u>)  N(w) 

Ht(w)'  H,(u>)X(w)’  Hi(w)lX{u>) 


Remark  2  The  expression  (B.2)  of  the  complex  cepstrum  reveals  that  the  received 
signal  is  not  minimum  phase  (the  complex  cepstrum  has  negative  components  in¬ 
troduced  by  the  noise).  This  illustrates  that  we  must  be  very  careful  when  we  want 
to  apply  inverse  filtering. 
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B.2.  DERIVATION  OF  THE  COMPLEX  CEPSTRUM  OF  THE 
MANNING- WINDOWED  CW  PULSE 

Let  us  now  c&lculate  the  complex  cepstrum  cf  the  transmitted  signal  *(<).  Here  x(t) 
is  a  Hanning- windowed  CW  pulse  given  by 

x(<)  =  0.5(1  -  cos  2xt/r)cos  wot, 

where  w0  is  the  CW  pulsation  and  T  is  the  length  of  the  Hanning  window.  The 
spectrum  X{w)  of  r(()  is 

X(u>)  =  (0.5<?o(w)  +  0.25Q0  (w  +  2t/2 T)  0.25(?o  (w  -  2jr/2T)] 

*  ^  ($ (w  —  wq  )  -f  6  (w  +  u'o )] , 


where 

^  ^  sin  uiT 

Q°(cw)  =  T— - — . 

Alternatively 

-Y(w)  =  0.25(?o(w  -  w0)  +  0.125Qo  -  w0  +  2x/2T) 

+  0.125Qo  (w  -■  w0  -  2ir/2T)  +  0.25<?0  (u>  +  w0) 

+  0.125(?0  (w  w0  4-  2x/2T)  +  0.125<?o  (u>  +  <*>0  -  2rj2T) . 


The  first  sidelobe  of  X(w)  is  quite  low  compared  to  the  principal  lobe  (X (first  lobe)/ 
•Y(w0)  =  0.00843)  and  thus,  one  does  not  make  a  serious  error  if  one  derives  the 
complex  logarithm  of  .Y(w)  from  only  the  principal  lobes  0  5Q0(w  -  w0)T  and 
0.5Qo(w  +  w0)7\  Under  this  assumption  X(w)  is  given  by 


*(«)  =  logf 


j,sin(w  -  w0)T  ^sin(w  +  w0)T 
(w  —  Wo  )T  (w  -f-  wo  )T 


We  now  apply  the  band-pass  mapping  defined  in  Subsubsect.  3.2.6.  We  recall  that 
it  consists  cf  band-pass  filtering  the  spectrum  X(w)  followed  by  the  mapping  trans¬ 
formation.  Here  the  spectra  is  filtered  around  the  frequencies  — wu  and  wo  in  such  a 
way  that 


-Yflit«red  (w 


*(w), 


if  0  <  |(w  -  wo)T|  <  ir 
otherwise. 


and  0  <  |(w  4-  wo)T|  <  t; 


■  B.2.I.  Derivation  of  X(u>)  for  0  <  |(w  -  w0)r|  <  w 

For  the  values  of  w  close  to  wo,  one  can  assume  that  the  quantity 


rsin(w  +  w0)r 
(w  4-  w 0)T 
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is  negligible  and  write 


(This  condition  is  fulfilled  when  w0  is  much  larger  than  l/T.)  If  one  now  applies 
band-pass  mapping,  the  frequency  u>  is  transformed  into  w  by 

<  _  w  -  (u>0  --  t/T) _ w  -  (w0  -•  t/T) 

“  "%  +  t/T  -  (w0  -  r/T)  ~  *  2t /T 


and  therefore 


w  ~z  2w  {T  -f  wo  —  t/T. 


The  expression  of  Ato0(w)  becomes,  with  respect  to  w  , 

v  /  \  t^sin(2w' -») 

*”»<" ( -  =T—j~nr' 

And  then  the  logarithm  can  be  expanded  in  series  as  follows: 

.  00  /  1\rt«>2n-l  □ 

X  ,(- )  =  log  ir  +  E  -  „(2-n)r  "(V  - 

where  Bn  are  the  bemouilii  numbers. 

k  B.2.2.  Derivation  of  A(w)  for  0  <  |(w  +w0)T\  <  t 

For  the  values  of  u>  close  to  -  wo,  one  can  assume  that  the  quantity 


,  sin(u>  -  u»0)r 

(w  —  u>o  )T 


is  negligible  and  write 


*—(«)  =  m*  ir  +  k,  . 


If  one  now  applies  band-pass  mapping,  the  frequency  w  is  transformed  into  u-  by 

'  _  w  —  (—Wo  +  w/T)  w  +  w*o  —  t/T 

W  r -w,,  -l-  t/T  -  (-u/o  -  w/T)  2 t/T 


and  therefore 


w  =  2w  /T  —  u’o  -f  w /T • 
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The  expression  of  X ^(w)  becomes,  with  respect  to  u> 


*-,(«•/)  =  £r 


sin(  -2u>  t) 

-2w'  T* 


Then  the  logarithm  can  be  expanded  into  series  as  follows: 


r_wo(u/,)  =  log  \T+Y, 

ns  1 


(-l)r,2»n-lBn 

n(2n)! 


(2w 


*)2n. 


The  sum  )  is  called  Xbp(w  ). 


■  B  2.3.  Inverse  Courier  transform  of  A’ep(w  ) 

Let  xpB(t)  be  the  inverse  fourier  transform  of  Xqp(u>  ),  given  by 

•^Bp(t)  =  J  Xbp(w  )eJ“  f  du) 

/o  ,  »»  , 

jf  _W0(w  )e;u'  ‘  du;  and  the  integral  /  Xu,0(o/)e,u' f  du/. 

/_„„  is  given  by 


J-wo  =  lo8 


*r  £  +  £  !=!£* ^  />'  +  *,\ 


And  if  the  variable  w  is  changed  into  the  variable  u;  =  2u>  +  the  integral  /_Wo 
becomes 

=2«->*'  log  IT  j  e’“"‘  d„ 

S"  /  w*"ejWl' 


+  E 

ns  1 


n(2n)! 


/ 


/«o 19  g»ven  ty 


-  “*  *7  [  +  E  ("lX)i'*‘  0 (2“'  - 
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And  if  the  variable  w  is  changed  into  the  variable  w  =  2u>  -2k  the  integral 
becomes  r 

log  \T  J  ejul/i  dw 

(-i)n23n— r 


/-a  =  2e^‘ 


-  (_i)"2a"->Bn  r 
+  nl  n(2n)! 


Consequently  zpn(t)  assumes  the  form 
ipB(0  =  4cos(^irt) 


logr.8^  +  £  r  ^v-/‘dw  . 

\*t  ^  n(2n)! 


(B3) 

This  can  be  simplified  as  follows.  When  I  =  I  wJnc>u,t/J 

dt  is  integrated  by  parts 

it  has  the  following  form: 

,  _  o J»t/i  .  4nfa;'n~*  _  8n(2n  -  l)w3n-3 

1  l  jt  +  t2  jt3 

(-l)n'l(2n)!2!n'V  (-  l)n-l(2r»)!22,,u>  (-l)n2ln+,i 

+  tln  +  J’ 

and  consequently  the  integral  /»  =  J  wineiut*1  dt  is  given  by 

_Jn+1  sin(f’irt)  4nir2"  cos(|-irf)  8n(2n  -  l)*-2"-1  sin(^-Tt) 

/.  =  IT  -j—  + 


t  4-ir< 
{2n- J 


t2 

4 

(-l)"-2(2n)!22n_1ir3  sin(frt)  (-l)"-‘(2n)!22,hr  cos(birt) 

1  + 


tJn-l 


\Wt 


or 


,  >.  8n(?n  -  l)w2n“l  (-l)r’-,(2n)!2:n-1>r2\ sin(^t) 

f'  =  'K' - o - - Jt 


+  ••  •  + 


ti  tin-l 

(-l)^-1(2n)!22^r  \  cos(^Tt)] 


nt 


f2n-l 


•  \  cos( 

V~F 


TkT 


And  if  we  define  the  variables  Cn(t)  and  Dn(t)  as 


Jn  8n(2n  —  l)*’2"'1  ( -l)"-2(2n)!22"-1ir2 


C„(t)  =  k 

l" 

„  .  „  4mr2n-1  {-l)"-‘(2n)!22nir 

^(0=— —  +  --+"Li>n-If - • 


<>«- J 
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The  simplified  form  of  £Pb(0  becomes 


( -l)n2,n_1  5n  ,  sin(isrt)  cos^xt)' 

iP9(0  =  4co.(^)-E  ~^)! - (C"{t)~]^T~  + 


\  cos(t<)  log  iTl 


B.S.  DERIVATION  OF  THE  COMPLEX  CEPSTRUM  OF  THE  BOTTOM  IMPULSE 
RESPONSE 

We  consider  Ht(u>)  as  the  transfer  function  of  one  finite-thickness  layer  given  by  (25) 


Ht(w)  =  -e 


(1  -  ^eJ"r)(l  +  y^e^r) 

(1  -  v^Te-i-Kl  -I-  v^e->“T)’ 


where  r  is  the  time  delay  across  the  layer  and  r0(r0  <  1)  is  the  reflection  coefficient 
of  the  first  layer  boundary.  One  must  normalize  ff;( <*■•)  by  removing  the  component 
,_e-jvr  W€  en(j  Up  w;th  the  simplified  form  of  Hi(uj): 

a- v^^t) 

U  '  (1  y/rSe-*")’ 

If  one  derives  the  complex  cepstrum  in  the  full  frequency  band,  the  following  ex¬ 
pression  is  obtained  for  hi(i): 

MO  =  -  r06(t  +  2r)  -  +  4r)  -  +  fir) 

+  •  •  +  r0#(<  -  2r)  +  *r^(f  -  4r)  +  -  6r)  +  •  •  • 

Under  the  assumption  that  the  terms  higher  or  equal  to  can  be  neglected,  one 
arrives  at  the  complex  cepstrum  expression  of  the  second  order: 

MO  =  -r06(t  +  2r)  -  +  4r)  +  r06(t  -  2r)  +  ^r^(t  -  4r). 


The  problem  is  that  the  bottom  transfer  function  has  been  band-pass  filtered  as  the 
Hanning- windowed  CW  pulse.  Consequently  the  same  band-pass  mapping  must  be 
used  before  deriving  the  complex  cepstrum.  We  recall  that  the  frequency  transform 
is  defin'd  by 

forO<|(u,-u.„)T|<*, 

|  x— (-^-/+y-/r)  for  0  <  |(u»  +  u/0)T|  <  x. 
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Let  us  now  derive  the  expression  of  Ht(u  )  around  the  two  frequencies  u>o  and  -u»#. 

■  B.3.1.  Derivation  of  for  0  <  |(ui  -  wo)T|  <  x 

Hi(w')  is  given  by 

_  [l  -  'T+u*-w/T)r]  [l  +  It+'*q-*IT)t) 

H'(U>  "  (1  -  v/fJ(-i(>“,/r+w»-./r)r]  jj  +  /»]. 

Let  us  call  a  the  quantity  w®  -  x/T,  and  take  the  logarithm  of  both  sides  of  the 
equation  and  expand  it  into  its  Taylor  series 

Ht{J)  =  -  s/Tlt**" /T*a)T  -  LroejI(I“,/r+e,)T 

-  ^Ov^eiS(,w  /T+a)T  ~  +  vV(2 J IT  +  o)r 

-  £r0e*,“'/T+a>T  +  |r0v^e>3^'/r+a>' 

+  •  •  •  +  'T+a*  -J-  lroe-^^’/r+ajr 

+  irov^e->3<,J/r+a)r - r+a)r 

+  ir0e~'',<Jw'/:r+a>r  -  ir0v/^e--»S(2tt',/r+a)r  +  •  •  • 

All  the  odd  terms  disappear  and  the  expression  becomes 

H,( u»)  =  -  roe>,<Jw'/T+°)T  -  lr2e^,1-'/r+a>r 

-  J.p3e>«(lw'/r+“)T  +  . . .  +  r0e~ji{iu' ,T+ a)r 

+  *r*«‘*(,w’/r+0,)r  +  ^Je"i#(,w,/r+o)r  +  • 


■  D.3.2.  Derivation  of  for  0  <  |{w  +  w0)T|  <  x 

Bv  the  same  kind  of  derivatioi  as  in  the  previous  case,  one  ends  up  with  the  following 
expression  for  ): 

H,{J)  -r0<.>*(*“'/T-a)r  _  ^P*c>«WI'-a>T 

-  Lp3eJ#<5“'/r-a>r  +  •  •  •  +  roe-^,“'/r  a)T 
+  fr2e-^,w'/r-a)r  +  /t—)t  x  . . . 
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■  B.S.3.  Derivation  of  the  complex  cepstmm  Aj(t) 
TcJcj  the  inverse  fouxier  transform  of  Hiiu  ) 


M«)  =  j  Hl(J)ei“tdJ 

“  j  ±  j  Ht(u>')eiu  tdu' 


and  derive  the  integrals  corresponding  to  the  terms  ro  of  since  the  derivation 

of  the  higher  order  terms(r,,  rj,  . . .  )  is  similar.  Define  first  the  function  j(t)  as 


/T-<n)rejJt  <JU>* 
/T-a)rejJt 
/Tra)r  ejJ  t  fa’ 
e~ji{iw  /T+a)re ju  t  fo' t 


which  gives 


g(t)  =  —-L—efl”  fl  -  e~J<‘-*'.'T) j 
jr(t-4r/T)  l  J 

+ - I _ e-’iat  LX'-W*)  -  ll 

+  j*(t  -  4r/T)e  le  lJ 

,€i*ar  j€i(‘+<f/J‘)  _  jj 


jw(t  +  4t/T) 
2* 

jx(t  +  4r/T)t 


-j  Jar 


[l  - 


r/T'} , 


or 


g(t)  =  2* 


cos(2qt) 


sinir(<  -  4t/T) 
w(t  -  4 t/~T) 


+ 


sin(2&r) 


sin*  ^»(t  -  4 r/T) 
~\x(t-4r/T) 


-  2r 


cos(2a7-) 


cosir(t  +  4  t/T) 
n (t  -r  4 t/T) 


sin(2ar) 


sin*  +  4t/T) 

H«  +  4r/D 
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Now  define  the  function*  hi{t  +  4 r/T)  and  ht(t  -  4 r/T)  as  follows: 


h,  (t  +  4 r/T) 


-2* 


cos(2ar) 


cos  r  (t  +  4t/T) 
*(i  4  4  r/T) 


sin(2ar) 


sin*  ^>r(t  4  4 t/T) 
^(t  +  4r/T) 


hj  (t  -  At /T)  =  2x 


cos(2or) 


sin  x  (t  -  4 r/T) 
w(t  -  At  IT) 


4  sin(2err) 


sin1  J-x  (t  -  4 t/T) 
\r(t-4r/T) 


The  derivations  of  the  integrals  of  higher  order  are  completed  in  the  same  way  and 
finally  the  complex  cepstrum  assumes  the  form 

MO  =  ro'hx  (t  4  4 r/T)  +  \r */»,  (<  4-  8 t/T)  +  (t  4  12r/T) 

4  •  •  •  4  r0ht  (t  -  4r/r)  +  fr’/i,  (t  -  8  r/T)  4  ^h,  (t  -  12  r/T)  +  ■■■. 

In  the  second  order,  the  complex  cepstrum  will  have  the  form 

MO  —  r0/»i  (<  4  4r /T)  +  (t  4  B~/T)  4  r0hj  (t  -  4 r/T)  4  (t  -  8 r/T) . 


Let  us  now  discuss  the  expression  for  the  complex  cepstrum  with  respect  to  the 
different  parameters  involved.  The  most  important  is  the  time  delay  r  whose  value 
determines  whether  or  not  we  will  be  able  to  separate  the  bottom  complex  cepstrum 
from  the  ‘boundary  reflections’,  r  is  obviously  a  function  of  the  layer  depth  and 
the  layer  sound  velocity  (i.e.  it  depends  on  bottom  composition).  Figure  112a, 
(18),  represents  the  sound  velocity  in  the  layer  function  of  the  porosity  of  the  layer 
components.  Figure  112b,  (19),  shows  the  dependence  of  the  sound  velocity  function 
on  relative  density,  porosity  and  reflection  coefficient.  Let  us  consider  two  relatively 
opposite  situations.  The  first  one  consists  of  a  layer  of  low  density,  high  porosity  and 
low  coefficient  of  reflection.  Let  us  take  the  case  where  these  three  parameters  have 
respectively  the  values  1.3,  80(%)  and  0.1.  According  to  Fig.  112b  the  corresponding 
velocity  equals  1400  m/s  if  one  assumes  a  water  sound  velocity  of  1500  m/s.  Let 
the  layer  depths  be  respectively  equal  to  50  in  and  100  m.  The  incidence  angle  0,b 
being  60°,  the  time  delays  r  are  respectively  0.041  s  and  0.082  s.The  corresponding 
complex  cepstra  M0  are  depicted  in  Figs  113a  and  b.  The  second  case  corresponds 
to  a  layer  of  high  density,  low  porosity  and  high  reflection  coefficient  with  parameter 
values  respectively  of  2.1,  32(%)  and  0.4.  Thus  the  sound  velocity  within  the  layer 
is  1800  m/s,  and  for  the  same  layer  thicknesses  and  incidence  angle  as  previously, 
the  time  delays  are  0.032  s  and  0.064  s.  The  complex  cepstra  M0  are  presented  in 
Figs.  114a  and  b. 

It  is  meaningful  to  assume  that  the  time  delay  4r  is  much  smaller  thar  t*  and 
therefore  M0  can  be  band-pass-filtered. 
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B.4.  DERIVATION  OF  THE  SURFACE  IMPULSE  RESPONSE  VFTER  BAND-PASS 
MAPPING 


We  consider  the  case  of  low  frequencies  and  a  very  smooth  •urfa  e.  Under  these 
assumptions  the  surface  transfer  function  assumes  the  form  ( Fraunhofer  diffraction), 
[20], 

-e'*, 


where 

=  — costf.C  (R„t,). 

c 

9,  is  the  incidence  angle  of  the  wave  with  respect  to  the  normal  to  the  surface; 

*s  the  surface  profile  at  the  point  of  specular  reflection.  is  a 

random  function  of  the  sea  surface  elevation  (surface  roughness).  Our  purpose  is  to 
derive  the  impulse  response  h.(t )  after  the  band-pus  mapping  has  been  applied.  As 
before  let  us  distinguish  the  two  cases  0  <  |w  -  u>0)T|  <  *  and  0  <  \w  -I-  w0)T|  <  x. 


■  B.4.1.  Derivation  of  H,(w)  for  0  <  \w  -  woJT’l  <  x 
The  frequency  transform  is  defined  by 

u>  —  2w  /T  +  u>q  -■  x /T. 


Thus, 


H.(w  )  =  exp 


-*( 


+  “><*  -  f  I  COS 


.2  o  .  _ 

’  1 

4w  .  _ 

exp 

-J —  cos  0,{(R„l,) 
c 

exp 

cos  <»,(,(  R„t,) 

where  a  =  w0  -  x/T. 


m  B  A. 2  Derivation  of  H,(w)  for  0  <  |w  +  u>o)7’|  <  x 
The  frequency  transform  is  defined  by 


Thus, 


.  2  /  2w' 

x  \ 

.  .1 

exp 

-  J  cos  9,((R„t,) 

» 

.  2o  „  „„ 

Aw 

1 

exp 

;  —  cos  9t({R„tt) 

c  J 

exp 
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Take  the  inverse  fourier  transform  to  get  h,{t) 


i,(t)  =  j°  H. 


.(w  )e}u  ‘  du  + 


f  H‘ 


(u>)eju  ‘du> 


By  replacing  H,(u>  )  by  its  values  h,(t)  becomes 


h.(t)  =  exp 
+  exp 


j  — cos(0,  )<(£,,  t,) 


2a 


cos(0,  )<(*.,<.) 


f  «XP  J^'  (t-  dw 

exp  |jw'  du» 


i.e. 

h,(t)  =  exp 


j  — cos  (0,)£(R,,t.) 
c 


-  exp 


2a 


-j-— co8(0,K(ft>ti,) 


[  1  -  exp  [-j*  (t  -  (4co8^/cr)C(fi„t,))] 
j(t  -  (4 cos 9,/cT)C(R„t,)) 

1  -  exp  (j>(<  -  (4cosfl./crK(fl,,f.))) 


j(t  -  (4  cos 9,/cT)((R„t,)) 


Finally  we  end  up  with 


h.(t)  =  cos 


—  cos(0,K(fi,,t,) 
c 


^sin[ir(<  -  (4 cos 0,/cT)((R,,t,))) 
ic{t  -  (4co96,/cT)((R,,t,)) 


f  sin 


—  cos  (0.)«R.,t.) 


^sin;  (^r(<-  (4  cos 9,/cT)((R,,t,))\ 
k*{t  -  (4cos0M/cT)((R„t,)) 


Let  us  assume  that  the  random  function  £( R,,t ,)  is  described  by  the  roughness 
parameter,  which  is  the  random  wave  height.  If  h  =  hohr,  where  ho  is  the  basic 
wave  height  and  hr  is  a  random  number  following  a  normal  distribution  with  rero 
mean  and  unit  standard  dcviaticn:  and  the  incidence  angle  is  equal  to  60°,  the 
surface  impulse  response  shown  in  Fig.  115  is  obtained. 

The  problem  now  is  to  properly  filter  the  complex  cepstrum  in  order  to  deconvolve 
the  wavelet  and  the  medium  response( boundary  reflection).  If  the  time  delay  r,  is 
small  with  respect  to  the  signal  *(t)  length,  then  h,(t  -  r,)  -  *  h,(t  -  2r,) 

will  be  overlapped  by  £(t)(in  our  example  the  time  delay  r,  is  equal  to  0.133  s).  It 
would  thus  be  reasonable  to  filter  in  such  a  way  that 

iu,(t)  =  *(()  +  hi{t),  (B.5) 


-58- 


SACLANTCEN  3M-30S 


i.€. 


MO  =  *i(e  -  (r,  -  T»)J  -  ^MO  *  A», [<  -  2(r,  -  rj))  +  hs(t  +  n) 

-  \hs (t)  *  hs(t  +  2r,)  -  hA(t  -  r,).  (B.6) 


iw(0  is  the  band- pass-filtered  cepstrum  around  t  —  0  and  represents  an  estimation 
of  the  complex  cepstrum  of  the  wavelet;  jm(0  is  equal  to  ,i(t)  -  sm(0  and  represents 
cn  estimation  of  the  complex  cepstrum  of  the  so-called  medium  response  (boundary 
reflection). 


Remark  3  The  complex  cepstrum  component  corresponding  to  the  medium  re¬ 
sponse  is  more  sensitive  to  the  noise  than  the  wavelet  component:  in  our  case  the 
noise  efect  on  the  deconvolved  wavelet  is  null. 

Let  us  assume  that  the  complex  cepstrum  is  altered  in  such  a  way  that  i„,(0  and 
sm(0  are  given  by  Eqs.  (B.5)  and  (B.6).  We  now  apply  the  inverse  homomorphic 
transform  to  get  ths  deconvolved  wavelet  and  the  deconvolved  medium  response. 


B.5.  DERIVATION  OF  THE  DECONVOLVED  WAVELET 
Take  the  fourier  transform  of  both  sides  of  Eq.  (B.5) 

SwM  =  x(uj)  +  £,m. 

Take  the  complex  exponential  of  both  aides  of  the  previous  equation 

Sw{u>)  =  X(w)R,H- 

Take  the  inverse  fourier  transform  to  get  the  deconvolved  wavelet  given  by 

M0  =  *(0*M0  (B.7) 
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B.8.  DERIVATION  OF  THE  DECONVOLVED  MEDIUM  RESPONSE 

Before  proceding  with  the  derivation,  we  stress  that  im(t)  is  strongly  perturbed  by 
the  noise  component  n(t)  (Eq.  (B.6)).  Since  we  recall  that  im(t)  is  given  by 

im(t)  =  M*  -  (T*  ~  ri)l  +  M*  +  n)  -  *  M*  ~  2(t»  -  r/)) 

-  \hz{t)  *  h3(t  +  2n)  -  h4(t  -  t,). 


Take  the  fourier  transform 

Sm(w)  =  4-  H3(u)e>“T>  - 

-  ^ff2(w)Vjw2(T*-T,)  -  Ht( u<)e-jur>. 

H3( u')e-’“'r|  -  ^/f3(w)3eJU,2r|  can  be  jonsidered  as  the  Taylor  expansion  of  log[l  + 
^3(u’)e',"Tl]  to  the  second  order.  /f2(w)e--'w(T* -T|*  -  ^if2(‘*')2e-3"2^T* _T>^  can  be 
considered  as  the  Taylor  expansion  of  log[l  +  /fj(u/)e ~ —  n)]  to  the  second  order. 
Thus,  Sm(u>)  will  have  the  form 

SmM  =  log(l  +  H2y)e~^r‘-^ 

+  log[l  +  H3(u>)e’wT>}  -  H4(u>)t-ju,T‘ . 


Take  the  complex  exponential  of  both  sides  of  the  previous  equation  and  reintroduce 
the  linear  phase  component  e~iWT> 

5m(w)  =[1  +  +  fl3(u;)] 

exp(-ff4(u')e_''“T*i. 


Let  Hp(u))  be  the  perturbation  exp ( —  (u^ )e “ ]  and  take  the  inverse  fourier  trans¬ 

form  of  the  previous  equation  to  get  the  deconvolved  medium  response  before  resti¬ 
tution  of  the  linear  phase  component 

Jm(0  =  ($(')  +  -  (r»  -  ri))  *  [*(<  -  d)  +  MO)  *  MO  (b.8) 
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B.7.  APPLICATION  TO  PARTICULAR  CASES 

We  now  consider  the  reflected  paths  to  be  solely  reflections  at  the  bottom  and 
bottom-surface  without  any  modification  of  the  tranfer  function  of  the  transmitted 
signal  up  to  a  time  delay.  In  other  words  the  bottom  transfer  function  and 

the  surface  transfer  function  H,( u/)  verify 

H,(u>)  =  1 

H,( u»)  =  1  for  anyO  <  u>  <  it. 

We  also  assume  that  the  signal-to-noise  ratio  is  high  enough  to  allow  all  terms 
containing  N (ui) / X (u>)  to  be  neglected.  Under  these  assumptions  the  relations  (B.7) 
and  (B.8)  assume  the  simplified  form 

■Su»(0  =  z(<)  (B.9) 

im(<)  =  (£(t)  r  i\t  -  (T,  -  r()]  *  6{t  -  r,).  (B.10) 


B.H.  CONCLUSIONS 

The  derivation  of  the  complex  cepstrum  using  this  simple  propagation  model  evokes 
two  major  observations: 

(1)  The  received  signal  after  reverberation  is  not  minimum  phase  (Eq.  (B  2)). 

(2)  The  derivation  of  the  layer  complex  cepstrum  shows  that  it  is  difficult  in 
the  cepstral  domain  to  separate  by  rectangular  windowing  the  transmitted 
pulse  from  the  layer  response,  because  they  occupy  the  same  cepstral  space 
(Figs.  I13a,b  and  114a,b).  Therefore,  what  .ve  deconvolve  is  the  wavelet,  from 
the  boundary  reflections,  and  incidentally  obtain  some  information  about  the 
layer. 

The  location  of  the  source  and  the  array  (closer  to  the  surface  or  closer  to  the 
bottom)  allows  us  to  estimate  the  bottom  impulse  response  or  the  surface  impulse 
response. 


01  - 


SACLANTCEN  SM-203 


Appendix  C 

Estimation  of  the  rank  of  the  correlation 
matrix  Rtl 


C.).  TEST  OF  AIC  AKAIKE 


The  method  developed  in  this  appendix  has  been  already  used  in  passive  array 
processing  [29].  The  purpose  was  to  estimate  the  number  of  sources  from  the  cross- 
spectral  matrix  measured  at  the  output  of  the  array 


Let  X  —  ( A'i  ,  X2 , A'/y)  be  a  series  of  independent,  zero-mean,  gaussian  vector 
random  variables  of  order  P  and  variance  matrix 

Rzx  ~  ai  +  R„. 

Their  probability  density  is  given  by 


n(A7ft,  R„)  =  2x“'v/,/!(dct«tx)"v/2exp(-i.  £  X*R^Xn). 

n=l 

Let  us  define  the  likelihood  fund  ion  by 


where 


${X/a,  R„)  =  \np(X/a,  R„)  -  Pln2ir 

N 


=  ^  det  Rxx  +  ^2  X  J Rxl  I„ 

n  —  l 

=  In  det{a/  +  R„)  +  +  R„)~l  R„ 


R 


xx  — 


1  N 


T 
n  ' 


One  wants  to  estimate,  in  the  maximum  likehood  sense,  the  two  unknowns  (a,  R,,)} 
which  is  equivalent  to  minimizing  the  function  $(X/a, 
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■  C.1.1.  Minimization  of  the  likelihood  function  i(X /a,  R,,) 

Let  us  assume  that  the  matrix  Rt,  is  of  rank  Q.  R„  can  be  expanded  into  its 
eigenvalues  and  eigenvectors  decomposition 

Q 

R „  =  V  A.u.  -if. 

«=•! 

The  inverse  of  ctl  +  R„  is  g'ven  by 

(al  +  R„)~l  =  —  I  -~7TU'U^  > 

or  z — 4  a  +  A, 

»= . 

and  then 

!  <5  A.  * 

tr(a/  +  R„)  lR x-  -  -  tr£«  -  V  — , 
a  a  +  Aj 


det(a/+«,.)  =  ap-£?fj(a  +  A,). 

«=i 

Thus  tiie  function  $(X/a,  R,,)  will  have  the  form 

P  1  T  Q  \ 

*(X/at\i,ui)  =  (P-Q)lna+  Vln(a  i- A,)  +  ^-  tr Rxx  -  T —±~uf  RxxUi  . 

• — '  2a  4—'  a  -f  A,- 

t=i  i=i 


To  globally  the  function  $(X/a,  A<,  u<)  minimize  one  minimizes  it  for  each  of  the 
three  variables  vdth  other  two  fixed. 

C.I.la.  Minimization  of  $  with  reaped  to  i  €  [!,<?],  with  Aif  »  €  jl,Q],  and  a  fixed 


This  is  equivalent  to  maximizing  the  quantity 


rc 


One  recalls  that  a  quadratic  form  uRxxuT.  ||  u  ||<  1,  is  maximum  if  u  is  the 
eigenvector  of  Rxx  corresponding  to  the  largest  eigenvalue.  Because  we  have  the 
sura  of  Q  quadratic  forms,  the  maximum  is  reached  ior  the  Q  eigenvectors  iq  of 
.ft**  Thus  the  minimum  of  the  function  $  is  given  by 

P  .  Q  .  . 

‘f(Z/a,Ai,u,)-(P-Q)ina  +  T'ln(a  +  A1)+—  ttRxx  -  V  ,  (C.i) 

2a  L  f^n  +  A'J 
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where  o',-, «  £  [?,<?],  are  the  eigenvalues  of  Rmm  arranged  in  decreasing  order. 


C  1.1b.  Minimization  of  $  with  respect  to  A ,,  »  £  [1,  Qj,  with  u,,»  6  [1,  Q],  and  a  fixed 
Noting  that 

A,  _  _  a 
(a  4  A«)  (a  +  A*) 


tr£„  =  &i 


it  follows  that 


f(X/a,  A,,u,)  =  [P  -  Q))no  +  £ln(a  4  A<)  +  ^<  +  j  ^ 


Q 

l  V'  <7’ 


2a4r+> 


The  function  is  minimum  for  the  gradient  equal  to  iiero  and  the  hessian  positive.  It 
can  be  easily  verified  that  the  hessian  is  always  positive.  The  gradient  equals  zero 
when 

(ITTaj  “  (q  +  a,-)*  =0,  1-,-P| 


which  gives  us  the  solution 


A i  -  &i  -  a. 


The  minimum  of  the  function  $  assumes  the  form 


*(*/a,A;,uO-(P-Q)lna+£ln(*)+^  £  *< + 


,  ,  Q 


C.l.lc.  Minimization  of  $  with  respect  to  a,  u<,  t  6  [1,Q]  with  A<,  i  £  (l,Q]  fixed  As 


previously,  $  is  minimum  for 


P-Q  i  .  n 

— -  -  „?  L  '<  =  °- 

•=<5+l 


and  the  solution  is  given  by 


1  - 
(P~Q)  Zv  <r" 
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The  global  minimum  of  is  given  by 


*(X/6,\itu,)  =  (P-  Q)\n 


E  *. 


i  (r  -  c?)  v 

i=Q  +  l 


where 


=  <£(<?)  +  —  +  LndetiJ«, 


W)MF  -0)lnfr„4^  £  *, 


(7* -2)  ^  1  ^ 

^7i  =  <J  +  l  J  »  =  <?  +  ! 


+  j><T,+  o 


1=1 


(C.2) 


-  £  In*,- 


C.2.  AKAfKE  CRITERION,  APPLICATION  TO  THE  ESTIMATION  OF  THE  RANK 
OF  R„ 

The  Akaike  estimate  of  the  order  of  a  mode]  at  the  minimum  cf  the  function  is 
,  „  .  .  number  of  free  parameters 

/(«)  =  -*(q)  + - ; 

where  $(q)  is  the  maxivnum  likehood  function  of  the  model  at  the  order  q.  In  our 
CMC  $(9)  is  given  by  Eq.  (C.2).  Using  the  Choleski  decomposition,  we  have 

RJt  =  LLt, 


where  L  is  ^  lower-trianguJor  matrix  with  P  rows  aru  Q  columns.  Then  the  number 
of  free  parameters  is 

P  +  (P  -  l)  +  (P  -  2)  +  •  •  (P  -  Q  +  1)  =  QP  -  \Q 7  =  Q(P-\Q), 


and  the  function  f(q)  takes  the  following  form 


f(Q)  =  (P-Q)  In 


1 

(P-Q) 


p 


E 

i=<?+l 


P 

-  ^  In  *•  + 

>=<3+i 


N 


in  which  the  constant  term  has  been  removed. 
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Fig.  1.  Three-path  propagation  model  (source  and  array  close  to  the 
surface). 


MD.T, 


Fig-  2.  Three-path  propagation  model  (source  aid  array  close  to  the 
bottom). 
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Fig.  3.  Canonical  representation  of  a  homomorphic  system. 


x(i)=x,(t).x,(t) 


B 

+  + 

D 

+ 

d* 

y(») 

X(1)--X,(1)VX2 

(») 

y(») 


Fig.  4.  Canonical  representation  of  a  homomorphic  deconvolution  system. 


Fig.  5.  Characteristic  system  D,  of  a  homomorphic  deconvolution. 


Fig.  6.  Bandpass  complex  cepstrum  system  D . . 
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Fig.  7.  Example  of  band-pass  mapping. 
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Fig.  8 a.  Global  complex  cepstrum  deconvolution. 


Fig.  8b.  Deconvolution  procedure  by  combination  of  homomorhic  and  Wiener 
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SHIFTEO  ANO  STRETCHED  POWER  SPECTRUM 


NORMALIZED  FREQUENCY 
Fig.  9.  Spectrum  of  the  CW  pulse  after  band¬ 
pass  mapping. 


NORMALIZED  FREQUENCY 


Fig.  10.  Spectrum  of  the  medium  response  after 
band-pass  mapping. 


NORMALIZED  FREQUENCY 


Fig  11  Spectrum  of  the  received  signal  after 
band-pass  mapping. 
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Fig  13.  Phase  of  the  medium  response  before  unwrapping. 


(aVrad) 
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Fig •  1 7.  First  derivative  of  the  phase  of  the  medium  response. 


Fig.  18.  Second  derivative  of  the  phase  of  the  medium  response. 


-  09  - 


SaCLANTCEN  SM-203 
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Fig.  1 9.  First  derivative  of  the  phase  oi  the  received  signal. 
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Fig,  20.  Second  derivative  of  the  phase  of  the  received  signal. 
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NORMALIZED  FREQUENCY 

'ig.  21.  Phase  of  the  CW  pulse  aftei  unwrapping  (tefore  linear  phase 
emovai). 


Fig.  22.  Phase  of  the  medium  response  after  unwrapping  (before  linear 
phase  removal). 
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?ig  23. 


Phase  of  the  received  signal  after  unwrapping  (before  linear 


these  removal). 


Fig.  24.  Phase  of  the  CW  pulse  after  unwrapping  (after  linear  phase 
removal). 
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Fig.  25.  Phase  of  the  medium  response  after  unwrapping  (after  linear 
phase  removal). 


NORMALIZED  FREQUENCY 

Fig.  26.  Phase  of  the  received  signal  after  unwrapping  (after  linear  phase 
removal). 
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Fig.  27.  Phase  of  the  received  signal  (explosive)  after  unwrapping  (before 
linear  phase  removal). 
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NORMALIZED  FREQUENCY 

Fig.  28.  Phase  of  the  received  signal  (explosive)  after  unwrapping  (after 
linear  phase  removal). 
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EIGENVALUE  OF  THE  AUTOCORR.  MATRIX 


CW  PARAMETERS 
RECTANGLE  W 
HALF  CYCLE  SIN  W. 
HANNING  WINDOW 
HAMMING  WINDOW 
LENGTH  16  SAMPL 
NORM  FREQ  :  0.25 
MAT.  ORDER  :  256 
S/M  -  3  0  <16 
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I  f  y  t  t  <  v  t  f  1  T  i  i  f  r+1  -r-E»  f  f  »  p  1‘l^r. 
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Fig.  29.  Eigenvalues  of  the  autocorrelation  matrix  of  the  received  signal 
(transmitted  pulse-  16  time-samples  CW  pulse). 


LxJ 

O 

3 


0l 


EIGENVALUE  OF  THE  AUTOCORR.  MATRIX 
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Fig.  30.  Eigenvalues  of  the  autocorrelation  matrix  of  the  received  signal 
(transmitted  pulse:  64  time- samples  CW  pulse). 
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Fig.  31.  Akaike  function  (transmitted  pulse:  16  time-samples  CW  pulse). 
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Fig  32.  Afcaike  function  (transmitted  pulse:  64  times-samples  CW 
pulse). 


SACLANTCEN  SM-303 


ESTIMATE  OF  THE  NOISE  TO  SIGNAL  RATIO 


Fig  33.  Estimate  of  the  noise-to-signai  ratio  (transmitted  pulse:  16  time- 
samples  CW  pulse). 


ESTIMATE  OF  THE  NOISE  TO  SIGNAL  RATIO 


Fig.  34.  Estimate  of  the  noise-to-signal  ratio  (transmitted  pulse:  64  time- 
samples  CW  pulse). 


-  017  - 


Saclantcbn  sm-jos 


BOnOM  REFLECTED  ARRIVAL 


Fig.  35.  Configuration  of  the  active  sonar  backscattering  simulation. 
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Fig.  3o.  Transmitted  pulse  (Hanning- windowed  CW  pulse). 
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Fig.  37.  Cross-correlation  between  the  transmitted  pulse  and  the  decon¬ 
volved  wavelet  (Hanning- windowed  CW  pulse;  3  multiples;  r  =  80,  110, 
170;  SNRR  =  14  <B). 
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Fig.  38.  Received  signal  (Hanning-windowed  CW  pulse;  3  multiples; 
r  =  80,  110,  170;  SNRR  —  14  dB). 


Fig.  39.  Received  signal  power  spectrum  (Hanning-windowed  CW  pulse; 
3  multiples;  r  =  80,  110,  170;  SNRR  =  14  dB). 
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Fig.  41.  Deconvolved  medium  response  by  copstrum  (Henning-windowed  CW 
pulse;  3  multiples;  r  =  80,  110,  170;  SNRR  =  14  dB). 
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Fig.  42.  Deconvolved  medium  response  by  Wiener  filtering  (Hanning-windowed 
CW  pulse;  3  multiples;  r  =  80.  110,  170;  SNRR  =  14  dB). 
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Fig.  43.  Received  signal  (Hanning-windowed  CW  pulse;  3  multiples;  v  =  80, 
110,  170;  SNRR  =  9  dB). 
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Fig.  44.  Received  signal  power  spectrum  (Hanning- windowed  CW  pulse; 
3  multiples;  r  =  80,  110,  170;  SNRR  =  9  dB), 
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Fig.  45.  Deconvolved  wavelet  (Hanning-windowed  CW  pulse;  3  multiples; 
r  =  80,  110,  170;  SNRR  =  9  dE). 
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Fig.  46.  Deconvolved  medium  response  by  Wiener  filtering  (Hanning-windowed 
CW  pulse;  3  multiples;  r  =  80,  110,  170;  SNRR  =  9  dB). 
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Fjg.  47.  Received  signal  (Hanning-windowed  CW  pulse;  3  multiples;  r  =  80, 
110,  170;  SNRR  =  4  dB). 
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Fig.  48.  Received  signal  power  spectrum  (Hanning-windowed  CW  pulse; 
3  multiples;  r  =  80,  110,  17C;  SNRR  =  4  dB). 
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Fig.  49.  Deconvolved  wavelet  (Hanning- windowed  CW  pulse;  3  multiples; 
r  -  80,  110,  170;  SNRR  =  4  dB). 
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Fig.  50.  Deconvolved  medium  response  by  Wiener  filtering  (Hanning-windowed 
CW  pulse;  3  multiples;  r  =  80,  110,  170;  SNRR  =  4  dB). 


Fig.  51.  Received  signal  (Hanning-windowed  CW  pulse;  3  multiples;  r  =  80, 
110,  170;  SNRR  = -1  dB). 
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Fig.  52.  Received  signal  power  spectrum  (Hanning-windowed  CW  pulse; 
3  multiples;  r  =  80,  1 .1 0,  170;  SNRR  =  -1  dB). 
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Fig.  53.  Ce convolved  wavelet  (Hanning- windowed  CW  pulse;  3  multiples; 
r  =  80,  110,  170;  SNRR  =  -1  dB). 
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Fig  54  Received  signal  (Hanning- windowed  CW  pube;  3  multiples  r  =  60 
110,  170;  SNRR  =  -6  dB). 
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Fig  55.  Received  signal  power  spectrum  (Hanning-windowed  CW  pube; 
3  multiples;  r  =  80,  110,  170;  SNRR  =  -6  dB). 
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Fig.  56  Deconvolved  wavelet  (Hanning- windowed  CW  pulse;  3  multiples; 
r  =  80,  110,  170;  SNRR  =  -6  dB). 
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Fig.  57.  Received  signal  (Hanning-windowed  CW  pulse;  3  multiples;  r  =  80, 
110,  170;  SNRR  =  -11  dB). 
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Fig.  58.  Received  signal  power  spectrum  (Hanning- windowed  CW  pulse; 
3  multiples;  r  =  80,  110,  170;  SNRR  =  -11  dB). 
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Fig.  59.  Deconvolved  wavelet  (Hanning-windowed  CW  pulse;  3  multiples; 
r  =  80,  110,  170;  SNRR  =  -11  dB). 
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Fig.  60.  Cross-correlation  between  the  transmitted  pulse  and  the  decon¬ 
volved  wavelet  (Hanning-windowed  CW  pulse;  5  multiples;  r  =  80,  110, 
140,  170,  200;  SNRR  =  16  dB). 
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Fig.  61.  Received  signal  (Hanning-windowed  CW  pulse;  5  multiples;  r  —  80, 
110,  140,  170,  200;  SNRR  =  16  dB). 
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Fig.  62.  Received  signal  power  spectrum  (Hanning- windowed  CW  pulse; 
b  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  16  dB). 
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Fig.  63.  Deconvolved  wavelet  (Hanning- windowed  CW  pulse;  5  multiples; 
r  =  80,  110,  140,  170,  200;  SNRR  =  16  dB). 
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Fig.  64.  Deconvolved  medium  response  by  cepstrum  (Hanning-windowed  CW 
pulse;  5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  16  dB). 
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Fig.  65.  Deconvolved  medium  response  by  Wiener  filtering  (Hanning- windowed 
CW  pulse;  5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  16  dB). 
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Fig  66.  Received  signal  (Hanning-windowed  CW  pulse;  5  multiples;  r  =  80, 
110,  140,  170,  200;  SNRR  =  11  dB). 
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Fig.  67.  Received  signal  power  spectrum  (Hanning-windowed  CW  pulse; 
5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  11  dB). 
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Fig.  70.  Deconvolved  medium  response  by  Wiener  filtering  (Hanning-windowed 
CW  pulse,  5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  11  dB). 


Fig.  71.  Received  signal  (Hanning- windowed  CW  pulse;  5  multiples;  r  =  80, 
110,  140,  170,  200;  SNRR  =  6  dB), 
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Fig.  72.  Received  signal  power  spectrum  (Hanning-windowed  CW  pulse; 
5  multiples,  r  =  80,  110,  140,  170,  200;  SNRR  =  6  dB). 


Fig.  73.  Deconvolved  wavelet  (Hanning- windowed  CW  pulse;  5  multiples; 
r  =  80,  110,  140,  170,  200;  SNRR  =  6  dB). 
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Fig.  74.  Deconvolved  medium  response  by  cepstrum  (Hanning- windowed  CW 
pulse;  5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  6  dB). 


Fig  75.  Deconvolved  medium  response  by  Wiener  filtering  (Hanning-windowed 
CW  pulse;  5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  6  dB). 
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Fig.  78.  Received  signal  (Hanning-windowed  CW  pulse;  5  multiples;  r  =  80, 
110,  140,  170,  200;  SNRR  =  1  dB). 
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Fig.  77.  Received  signal  power  spectrum  (Hanning- windowed  CW  pulse; 

5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  1  dB). 


RECEIVED  SIGNAL  POWER  SPECTRUM 


-<939- 


SACLANTCEN  SM-303 


fig  78.  Deconvolved  wavelet  (Hanning-windowed  CW  pulse;  5  multiples; 
r  =  80,  110,  HO,  170,  20C;  SNRR  =  1  dB). 


Fig.  70.  Deconvolved  medium  response  by  cepstrum  (Hanning- windowed  CW 


poise;  5  multiples;  r  =  80,  110,  140,  170,  200;  SNRR  =  1  dB). 
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Fig.  80.  Deconvolved  medium  response  by  Wiener  filtering  (Hanning- windowed 
CW  pulse;  5  multiples;  r  =  80.  110,  140,  170,  200;  SNRR  =  1  dB). 
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Fig.  81.  Received  signal  (Hanning-windowed  CW  pulse;  5  multiples;  r  =  80, 
110,  140,  170,  200;  SNRR  =  -4  dB). 
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Fig.  82.  Recalved  signal  power  spectrum  (Hanning- windowed  CW  pulse; 
5  multiples;  r  =  80,  110,  HO,  170,  200;  SNRR  =  -4  dB). 
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Fit  83.  Deconvolved  wavelet  (Hanning-windowed  CW  pulse;  5  multiples; 
r  =  80,  110,  140,  170,  200;  SNRR  =  -4  dB). 
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Fig.  84.  Configuration  of  the  active  sonar  backscattering  experiment 
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85.  Transmitted  pulse:  Hanning-windowed  CW  pulse  of  2  s. 


Fig.  86.  Received  signal  (beam  6,  ping  9,  range  8). 
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Fig,  88.  Deconvolved  medium  response  by  Wiener  filtering  with  the  decon¬ 
volved  wavelet  as  input 
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Fig.  89.  Deconvolved  medium  response  by  Wiener  filtering  with  the  trans- 
mitted  pulse  as  input 


Fig.  90.  Configuration  of  the  explosive  experiment. 
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Fig.  91  Received  sign  el  on  each  hydrophone  of  the  vertical  array 
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Fig.  93.  Received  signal  on  hydrophone  17. 


Fig-  93.  Power  spectrum  of  the  received  signal  on  hydrophone  17. 
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Fig  94.  Received  signal  on  hydiophone  4. 
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Fig.  97.  Comparison  of  the  partial  energies  of  the  decon¬ 
volved  wavelet  on  hydrophone  17  and  deconvolved  wavelet 
on  hydrophone  4. 
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Fig.  98.  Power  spectrum  of  the  deconvolved  wavelet  on  hydrophone  17. 


Fig  99.  Power  spectrum  of  the  deconvolved  wavelet  on  hydrophone  4. 
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Fig.  100.  Deconvolved  medium  response  on  hydrophone  17  by  Wiener  filter 
ing. 
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Fig.  10J.  Deconvolved  medium  response  on  hydrophone  4  by  Wiener  filtering. 
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Fig,  102 a.  Deconvolved  medium  response  on  each  hydrophone 
of  the  vertical  array  by  Wiener  filtering. 
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Fig.  103.  Zero-pole  diagram  of  the  transfer  function  G,r{z). 
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Fig.  104.  Zero-pole  diagram  of  the  transfer 
function  A'o(r). 
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Fig.  105.  Zero-pole  diagram  of  the  transfer  function  Jf|(z). 
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Fig  109.  Mixed  phase  signal 
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Fig.  110.  Comparison  of  the  partial  energies  of  the  signals  z0(f).zt(<)  and 
*>(«)• 
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Fig.  111.  Expression  of  the  complex  cepstrum. 
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Fig  112 a.  Sound  speed  vi  porosity  in  underwater  sediments 
made  on  core  samples  (from  Urick  [18]). 
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Fig.  112b.  Relationship  between  the  main 
characteristics  of  underwater  sediments. 
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Fig.  113.  Layer  complex  cepetram  with  porosity  =  80%:  (a)  layer  depth  = 
50  m;  (b)  layer  depth  =  100  m. 
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Fig.  114.  Layer  complex  cepstrum  with  porosity  --  32%:  (a)  layer  depth  = 
50  m;  (b)  layer  depth  =  100  m. 
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Fig.  115.  Surface  impulse  response. 
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